"""
Part of comet/snmr/modelling
"""
# COMET - COupled Magnetic resonance Electrical resistivity Tomography
# Copyright (C) 2019 Nico Skibbe
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <https://www.gnu.org/licenses/>.
import numpy as np
import matplotlib.pyplot as plt
try: # pg.__version__ >= 1.1
from pygimli.viewer.mpl import drawModel1D
except ImportError: # pg.__version__ < 1.1
from pygimli.mplviewer import drawModel1D
[docs]def brooksCorey(z, water_table, porosity, lam=1.6,
height_zero=0.12):
"""
after Brooks and Corey (1964)
\cite{}
"""
saturation_eff = np.ones_like(z)
vadose = np.where(z >= water_table + height_zero)
height = z[vadose] - water_table
saturation_eff[vadose] = (height_zero/height)**lam
return saturation_eff
[docs]def archie(porosity, saturation, water_resistivity,
tortuosity=1.0, cementation=1.3,
saturation_exponent=2.0,
formation_factor=None):
"""
porosity(z), saturation(z)
returns resititvity_bulk(z)
\cite{}
"""
if formation_factor is None:
formation_factor = tortuosity / (porosity**cementation)
# print('Archie: formation factor: {:2.2f}'.format(formation_factor))
return water_resistivity * \
formation_factor / (saturation**saturation_exponent)
[docs]def effectiveSaturationToWater(saturation_eff,
water_saturation,
water_residual=0.05):
"""
saturation_eff = (water - water_residual)/
(water_saturated - water_residual)
"""
water = saturation_eff * (water_saturation - water_residual) +\
water_residual
return water
[docs]def costabel(saturation, t2_saturation, lam=1.6):
""" \cite{costabel2011NSG}
Costabel, S., and U. Yaramanci, 2011, Relative hydraulic conductivity
and effective saturation from Earth’s field nuclear magnetic
resonance – a method for assessing the vadose zone: Near Surface
Geophysics, 9, 155–167.
"""
return t2_saturation * saturation ** (1./lam)
[docs]def modelVadose(z,
water_table,
porosity,
t2_saturated,
water_resistivity,
height_zero=0.12,
water_residual=0.05,
lam=1.6,
verbose=False,
**kwargs):
"""
Calculates a synthetical vadose zone on basis of a Brooks-Corey model
for saturation over the vadose zone, whereas lambda is the pore size
distribution index.
Also calculates the electrical resistivity(z) via Archies law, as well as
the distribution of relaxation times based on Costabel and Yaramanci
(2011).
returns (z, resistivity, water_content, relaxation_times)
"""
sat_eff = brooksCorey(z, water_table, porosity, lam=lam,
height_zero=height_zero)
water_content = effectiveSaturationToWater(
sat_eff,
water_residual=water_residual,
water_saturation=porosity)
saturation = water_content / porosity
resistivity = archie(porosity, saturation, water_resistivity, **kwargs)
relaxation_times = costabel(saturation, t2_saturated, lam=lam)
if verbose:
print('resistivity [$[\Omega$m$]$] (min/max): {:2.2f}, {:2.2f}'
.format(np.min(resistivity), np.max(resistivity)))
print('saturation [1] (min/max): {:2.2f}, {:2.2f}'
.format(np.min(saturation), np.max(saturation)))
print('water content [1] (min/max): {:2.2f}, {:2.2f}'
.format(np.min(water_content), np.max(water_content)))
print('relaxation times [s] (min/max): {:2.2f}, {:2.2f}'
.format(np.min(relaxation_times), np.max(relaxation_times)))
fig, ax = plt.subplots(1, 4, sharey=True)
drawModel1D(ax[0], depths=z, values=saturation,
xlabel='saturation [1]')
drawModel1D(ax[1], depths=z, values=resistivity,
plot='semilogx')
drawModel1D(ax[2], depths=z, values=water_content,
xlabel='water content [1]')
drawModel1D(ax[3], depths=z, values=relaxation_times,
xlabel='relaxation times [s]', plot='semilogx')
ax[0].set_ylim(min(z), max(z))
return saturation, (resistivity, water_content, relaxation_times)
[docs]def test_local():
z = np.arange(0, -4., -0.1)
water_table = -3.0
porosity = 0.3
water_resistivity = 20.0
t2_saturated = 0.2
water_residual = 0.05
lam = 1.6
height_zero = 0.12
saturation, (resistivity, water_content, relaxation_times) = modelVadose(
z,
water_table,
porosity,
t2_saturated,
height_zero=height_zero,
water_residual=water_residual,
lam=lam,
water_resistivity=water_resistivity,
verbose=True)
if __name__ == '__main__':
test_local()
pass
# The End