ARTICLE IN PRESS
Computers & Geosciences 34 (2008) 1080–1092 www.elsevier.com/locate/cageo
A hydrogeophysical synthetic model generator$ Bernard Giroux, Michel Chouteau De´partement des ge´nies civil, ge´ologique et des mines, E´cole Polytechnique de Montre´al, C.P. 6079, Succ. Centre-Ville, Montre´al, Que., Canada H3C 3A7 Received 20 January 2007; received in revised form 6 July 2007; accepted 15 November 2007
Abstract HGmod is a computer program that builds on stochastic realizations of porosity fields to derive electrical conductivity, dielectric permittivity and hydraulic permeability models. The presence of clay, the influence of salinity as well as temperature of the fluid of imbibition are taken into account in the underlying formulations. The saturated and unsaturated zones are also considered through the application of a saturation profile on the porosity field. A microgeometrical model is used to relate the porosity to the clay fraction. This model also used to derive an expression for the pore specific surface of the sand–clay mixture. The specific surface is subsequently used to compute the conduction at the pore surface when building electrical conductivity models. Dielectric permittivity fields are built by successive applications of either the Hanai–Bruggeman or Maxwell–Garnett mixing models, depending on the relative proportions of sand, clay, water or air. In addition, the dielectric permittivity of water and clay follow a Cole–Cole behavior. HGmod is therefore a versatile tool useful to generate synthetic datasets needed to anticipate the geophysical response under specific conditions and to study hydrogeophysical sensitivity or resolution analysis. r 2008 Elsevier Ltd. All rights reserved. Keywords: Hydrogeophysics; Modeling; Porosity; Electrical conductivity; Dielectric permittivity
1. Introduction The shallow subsurface is an important reservoir for groundwater supply. This zone is also used for the storage of all types of wastes, and is also exposed to many kinds of contaminants. Groundwater flow and contaminant transport are much dependent on the hydrogeological properties of the medium that are usually complex and scale $
Source code available at http://geo.polymtl.ca/giroux/ HGmod/. Code available from server at http://www.iamg.org/ CGEditor/index.htm. Corresponding author. Tel.: +1 514 340 4711x5233; fax: +1 514 340 3970. E-mail address:
[email protected] (B. Giroux). 0098-3004/$ - see front matter r 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.cageo.2007.11.006
dependent. The discipline of hydrogeophysics has emerged in the recent years due to the growing need to characterize and monitor the shallow subsurface in a non-invasive manner. Particular fields of study currently addressed are the improvement of the spatial resolution admissible by the inversion of geophysical data and its connection with the scale of variation the underlying physical parameter (Day-Lewis et al., 2005; Hubbard et al., 1999), the sensitivity of the results to variations in the physical properties of the subsurface (Chang et al., 2006; Furman et al., 2004) and the integration of geophysical and hydrogeological measurements (Tronicke and Holliger, 2005).
ARTICLE IN PRESS B. Giroux, M. Chouteau / Computers & Geosciences 34 (2008) 1080–1092
It is now a common practice for researchers to test their hypothesis numerically before performing costly measurements in the field. The usual numerical earth model used by geophysicists is a construction of the spatial distribution of the physical parameter measured by the technique employed. For example, a dielectric permittivity field with given spatial properties is generated to model the propagation of radio waves, or a model of resistivity is built to simulate an electrical resistivity tomography (ERT) survey. This way of doing has one serious limitation in that the physical property of interest is only weakly related to the underlying hydrogeological parameters (porosity or permeability). It is therefore difficult to assess quantitatively the sensitivity of a given geophysical method to the variation of the hydrogeological parameters. In this paper, we present HGmod, a program that allows the generation of a 2D or 3D distribution of dielectric permittivity, electrical resistivity and hydraulic permeability, all derived from a stochastic realization of a porosity field. These models can be used with other modeling programs to generate synthetic datasets. The presence of clay and the influence of salinity and temperature of the fluid of imbibition are taken into account for the calculation of the dielectric permittivity and electrical resistivity. The saturated and unsaturated zones are considered through the application of a saturation profile on the porosity field. For instance, HGmod can serve to generate different models for one varying parameter (e.g. salinity or water table level), in order to generate as many synthetic datasets that can in turn be used for sensitivity or resolution analysis. 2. Methodology 2.1. Building a porosity model Most shallow aquifers are made of unconsolidated granular material with interstice permeability, as opposed to fractured rocks with fracture permeability. Unconsolidated sediments present stochastic spatial structures, and their porosity f is adequately represented by a normal distribution (Freeze, 1975; Sun and Koch, 1998). Geostatistical simulations allow one to generate multi-dimensional random fields. The FFT-MA simulation algorithm (Le Ravalec et al., 2000) was implemented to generate 2D or 3D
1081
porosity fields that will serve as the basis for the calculations. This algorithm was retained for its simplicity, efficiency and ease of implementation. It can handle any covariance as long as the field is stationary (the probability distribution at a fixed position is the same for all positions, a weak assumption in this case). This generator, as its name implies, results from the combination of the moving average method with the FFT algorithm. The FFTMA generator relies on randomness components drawn in spatial domain. The latter, together with the discrete covariance term, are efficiently transformed to the spectral domain with the FFT algorithm, then merged and inverse transformed to yield the simulated field. One appealing feature is that the structural parameters are separated from the random ones. Using the latter feature, different realizations can efficiently be produced using the same covariance term and perturbing only the random one. We rely on this feature when building more that one model. With the FFT-MA algorithm, however, it is not possible to directly condition the simulated field to known data. Nevertheless, this can be done a posteriori using simple kriging. A basic implementation of simple kriging is included in the code to allow this. Given S the non-conditional simulated porosity field and conditioning data Z, the conditional simulated field T is (Chile`s and Delfiner, 1999) T ¼ S þ ½Z S ,
(1)
where superscript stands for simple kriging estimator. Using this approach, the conditioning values are exactly recovered in the simulated field. The quality of the realizations obtained with our implementation of the FFT-MA algorithm is succinctly examined through one example of a non-conditional simulation. The covariance model is composed of an isotropic spherical model with range 50 and sill 0.9 plus a nugget effect with sill 0.1. The domain has a size of 1024 1 1024 (2D models are generated by setting a single grid node in one of the three spatial dimensions). One hundred realizations are drawn, for which the average statistics are given in Table 1. The means of means and variances are close to their theoretical values (0 and 1, respectively) and dispersion is small. The simulated univariate distributions are then compared to the expected standard Gaussian distribution via a quantile–quantile plot (Fig. 1a). Also, the bivariate distributions are examined through the simulated variograms (Fig. 1b). In each case, one
ARTICLE IN PRESS B. Giroux, M. Chouteau / Computers & Geosciences 34 (2008) 1080–1092
1082
significantly affect both hydrogeological and geophysical properties. At this point, it is interesting to see how clay content can be taken into account when building a hydrogeophysical model. Marion et al. (1992) proposed a micro-geometrical model for mixtures of sand and clay in which clay is dispersed in the pore space of the sand matrix. According to this model, as clay packets fill the pore space, the porosity can be related to the clay fraction Cl by
observes that the average simulated statistics match very well the theoretical model. Fig. 2 shows four realizations of the porosity field with four different covariance models. In all cases, the mean porosity is 0.25, and the same seed was used to generate the random term in the simulations. Exponential, hyperbolic and stable covariance functions with ranges ax ¼ 10 m and az ¼ 3 m and a sill C ¼ 4 104 were, respectively, used in Figs. 2a–c, whereas a Gaussian covariance function with the same ranges and sill and a small nugget effect (C ¼ 105 ) were used to produce Fig. 2d (see Chile`s and Delfiner, 1999 for covariance function definitions). Note the similarities in the overall trends of these realizations. This is due to the decoupling of the structural parameters from the random ones in the FFT-MA algorithm and the use of the same random term for all simulations. It should be mentioned that although only 2D models are presented here, 3D models can be built with our program as well.
f ¼ fs Clð1 fc Þ,
where fs is the porosity of pure sand and fc is the porosity of pure clay. The clay fraction Cl is defined here as the volume of room dry clay (i.e. with associated bound water and macro-porosity) over the volume of room dry sand–clay mixture (Marion et al., 1992). By specifying values of fs and fc , it is therefore possible to compute a value of Cl at every porosity point, such as shown in Fig. 3 (fs ¼ 0:35 and fc ¼ 0:45 were used here). Note that Eq. (2) is valid as long as Clofs . We have restricted the domain of application of our program to this limit because soils with larger clay fraction are generally of limited hydrogeological interests due to their low permeability.
2.1.1. Including a clay fraction It is rather common that granular aquifers include a fraction of silt or clay. Clay or silt Table 1 Average statistics of 100 isotropic realizations Mean of means
Variance of means
Mean of variances
Variance of variances
0.001235
0.001147
1.011419
0.001236
(2)
2.2. Electrical conductivity model The presence of clay implies that the widely used formula of Archie (1942) cannot be used to model
QQ plot of Simulated Data vs Std Normal
Variogram of Simulated Data 1.2
2.0 1.0
1.0
0.8
0.5
Variogram
Simulated Quantiles
1.5
0.0 -0.5
0.6 0.4
-1.0 -1.5
0.2
-2.0 0.0 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 Standard Normal Quantiles
0
10
20
30
40 50 60 Distance [m]
70
80
90 100
Fig. 1. Validation of reproduction of univariate and bivariate distributions, from 100 non-conditional realizations.
ARTICLE IN PRESS B. Giroux, M. Chouteau / Computers & Geosciences 34 (2008) 1080–1092
0
1083
0
0.31
0.31
0.29
0.29
5
5 0.27
0.27 Z [m]
0.25 0.23
15
10
0.25
Z [m]
10
0.23
15
0.21
0.21
20
20 0.19
25
0.19 25
0.17 0
5
10
15
20
25
0.17 0
5
10
X [m]
15
20
25
X [m]
0
0
0.31
0.31
0.29
0.29 5
5
0.27
0.27 Z [m]
0.25 0.23
15
10
0.25
Z [m]
10
0.23
15
0.21
0.21 20
20
0.19
0.19 25
0.17 0
5
10
15
20
25
25
0.17 0
5
10
X [m]
15
20
25
X [m]
Fig. 2. Four realizations of porosity field: (a) Experiental covariance, (b) hyperbolic covariance, (c) stable covariance with a ¼ 1=2 and (d) Gaussian covariance with nugget effect. 0
0.32 0.28
5
0.24 10 Z [m]
the electrical conductivity of the medium. This is due to the increased surface conduction that occurs in media with high specific surface such as clay. Many models were developed to take into account the effect of surface conduction (Waxman and Smits, 1968; de Lima and Sharma, 1990; Pride, 1994; Revil, 1999; Tabbagh et al., 2002). In this work, the macroscopic model of Pride (1994) is used for that purpose. Pride derived his equations by using averaging of the equations known to apply at the fluid and solid phases, i.e. the microscopic scale, over a volume larger than the solid grains. The underlying assumptions are: the fluid is an ideal electrolyte (concentration less than 1 mol/l), isotropic model, and the thickness of the doublelayer is much less than the radii of curvature of the solid grains. According to this model, the bulk
0.20 15 0.16 20
0.12 0.08
25 0
5
10 15 X [m]
20
25
Fig. 3. Clay fraction model derived from porosity field of Fig. 2a.
ARTICLE IN PRESS B. Giroux, M. Chouteau / Computers & Geosciences 34 (2008) 1080–1092
1084
conductivity s is given by (Pride, 1994; Carcione et al., 2003) f 2½C em þ ReðC os ðoÞÞ n sf S þ ðS=mÞ, (3) s¼ T L where f is the porosity; T the tortuosity (for negligible surface conductivity, the ratio f=T is equal to the inverse of the electrical formation factor (Pride et al., 1993)); sf the fluid (ionic) conductivity; S the saturation; n the saturation exponent; C em the excess conductance associated with the electro-migration of double-layer ions; C os the frequency dependent conductance due to electrically induced streaming of the excess double-layer ions (coined ‘‘electro-osmotic’’ by Pride); L a geometrical parameter roughly equal to twice the ratio of volume of pores over surface of pores (Johnson et al., 1986; Johnson and Sen, 1988). We have introduced the saturation S in the original equation of Pride (1994) by multiplying sf with the saturation index of the second Archie equation (Scho¨n, 2004). Our formulation is valid as long as the saturation is high enough to maintain the double-layer, which is a reasonable working hypothesis due to the extremely thin thickness of the double-layer. The conductivity of a fluid with L ionic species is (Pride, 1994) sf ¼
L X ðezl Þ2 bl Nl
ðS=mÞ,
S¼
Sp Vh
ðm1 Þ,
(7)
where Sp is the surface of the pores and V h is the volume of the host matrix. From the definition of porosity (f ¼ V p =V with V p the volume of pores and V the total volume), Eqs. (2) and (7), it can be shown that S ¼ Cl Sc þ ð1 ClÞSs ,
(8)
where Sc is the specific surface of pure clay and Ss is the specific surface of pure sand. Then, from the definitions of specific surface and porosity and approximating L 2V p =S p , we get L
2f Sð1 fÞ
ðmÞ.
(9)
(4)
l¼1
2.3. Dielectric permittivity model
where zl is the valence (ezl represents the net charge and sign of the ion, e being the electric charge 1:6 1019 C), bl is the mobility and Nl is the bulk ionic concentration. The mobility can be approximated by 1=ð6pZf Rl Þ with Zf the fluid viscosity and Rl the effective ion radius. The viscosity of water is strongly dependent on temperature. For temperature in the range of 0–20 1C, the viscosity in centiPoise is (Weast and Astle, 1980) Zf ¼ 100 exp 2:303
Note that the effect of salinity on viscosity is neglected. Values of effective ion radius R for sodium and chlorine are, respectively, 0.116 and 0.167 nm. The bulk ionic concentration Nl (in m3 ) is equal to 1000sN a with s the salinity in mol/l and N a Avogadro’s number ð6:022 1023 mol1 Þ. Details on the computation of C em and C os are given in Appendix. To derive parameter L, we extend the model of Marion et al. (1992) to take into account the specific surface, itself related to L. The specific surface S is defined here as
Druchinin (2000) proposed models to compute the dielectric permittivity of sandy–clayey soils with various clay contents. The soil is considered to consist of four phases: air, water, clay and sand, in volume proportions of fa , fw , fc and fs , respectively. Depending on the proportion of each constituent, the models consist in successive applications of Hanai–Bruggeman or Maxwell–Garnett
1301 3:30233 998:333 þ 8:1855ðT 20Þ þ 0:00585ðT 20Þ2
and for the range 20–100 1C Zf ¼ 1:002 exp 1:3272ð20 TÞ 0:001053ðT 20Þ2 2:303 . T þ 105 (6)
(5)
models. In the mixing model of Hanai–Bruggeman, the permittivity of the mixture m can be obtained from (Scho¨n, 2004) W i m h 1 fi ¼ , i h m
(10)
ARTICLE IN PRESS B. Giroux, M. Chouteau / Computers & Geosciences 34 (2008) 1080–1092
where h is the permittivity of the host medium, i is the permittivity of the inclusion, fi is the volume fraction of the inclusion and W ¼ 13 for spherical inclusion. The term of interest m in Eq. (10) cannot be isolated. To solve the system, we use a modified version of Powell’s method (Press et al., 1992) to minimize the norm W i m h . 1 fi i h m In the Maxwell–Garnett model, the permittivity of the mixture is (Sihvola, 2000) fi P i h =h þ N j ði h Þ 3 j¼x;y;z , m ¼ h þ h P f 1 i j¼x;y;z N j ði h Þ=h þ N j ði h Þ 3
(11) where N j is the depolarization factor in j direction. The depolarization factors for flattened ellipsoidal inclusions with axes ð1; x; xÞ are used in our specific case (for details see Sihvola, 2000). The following four scenarios proposed by Druchinin (2000) are considered to compute the effective permittivity of the soil. In the first two cases, there is no clay in the model. The first scenario corresponds to a low moisture situation (fw o0:11). In this case, the procedure is as follows: (1) use the Maxwell–Garnett model to compute the permittivity of the sand and water
1085
mixture with water taken as the inclusion material; (2) use the Hanai–Bruggeman model to compute the permittivity of mixture 1 and air where air is considered the inclusion material. The second procedure, for (fw 40:11), is the same as above except that the Hanai–Bruggeman model is used in the first step. In the last two cases, the amount of clay is considered small. In the following, mixture 1 stands for the water/clay mixture, and mixture 2 stands for the combination of mixture 1 with sand. The third scenario is when fc þ fs ofw , and the following steps are applied: (1) use the Maxwell–Garnett model to compute the permittivity of the water and clay mixture with clay taken as the inclusion material; (2) use the Maxwell–Garnett model to compute the permittivity of the sand and mixture 1 combination with mixture 1 taken as the inclusion material; (3) use the Hanai–Bruggeman model to compute the permittivity of mixture 2 and air where air is considered the inclusion material. The last case is when fc þ fs 4fw , and the steps are as following: (1) use the Hanai–Bruggeman model to compute the permittivity of the clay and water mixture with water taken as the inclusion material;
Permeability [mD]
Permeability [mD]
0
175 150
5
70
0
60
5
50
125 10 100 15
Z [m]
Z [m]
10 40 15 75
20
50 25
25 0
5
10
15
20
25
30 20
20 10
25 0
5
10
15 X [m]
X [m]
Fig. 4. Hydraulic permeability fields derived from porosity model of Fig. 2a.
20
25
ARTICLE IN PRESS B. Giroux, M. Chouteau / Computers & Geosciences 34 (2008) 1080–1092
1086
σ [S/m]
σ [S/m] 0.020
0
0.031
0
0.019 5
0.030 5
0.018
0.029
0.017 Z [m]
0.016 0.015
15
10
0.028
Z [m]
10
0.027
15
0.014
0.026
0.013
20
20 0.025
0.012 0.011
25 0
5
10
15 X [m]
20
0.024
25
25
0
5
Re(ε)
10 15 X [m]
20
25
Re(ε) 11.0
0
19.0
0
18.5
10.5 5
17.5
9.5 9.0
15
17.0
10
16.5
Z [m]
Z [m]
10
18.0
5
10.0
16.0 15
8.5
15.5 15.0
8.0 20
20
14.5
7.5
14.0
7.0
25 0
5
10 15 X [m]
20
13.5
25
25
0
5
Im(ε)
10 15 X [m]
20
25
Im(ε)
0
0.90
0
5
0.75
5
10
0.60
10
5.90 5.85
Z [m]
Z [m]
5.80
0.45
15
5.75 15 5.70
20
25 0
5
10 15 X [m]
20
25
0.30
20
0.15
25
5.65 5.60 0
5
10 15 X [m]
20
25
Fig. 5. Influence of clay content on electrical conductivity and relative dielectric permittivity: (a), (c), (e) Clean sand and (b), (d), (f) Sand with clay content.
ARTICLE IN PRESS B. Giroux, M. Chouteau / Computers & Geosciences 34 (2008) 1080–1092
(2) use the Hanai–Bruggeman model to compute the permittivity of mixture 1 and sand, where sand is considered as the inclusion material; (3) use the Hanai–Bruggeman model to compute the permittivity of mixture 1 and air, where air is considered as the inclusion material. Druchinin (2000) suggests using a value of 35 for the effective axial ratio x to compute the depolarization factor when applying the Maxwell–Garnett model. In our implementation, the effective dielectric permittivity of water w and of clay c follow a Cole–Cole behavior: 0 1 is þ , 1 ðiotÞq o
(12)
where 1 is the optical dielectric permittivity, 0 is the static dielectric permittivity, t is the relaxation time and q is the Cole–Cole distribution parameter (Hasted, 1973). Using this model, dispersion is considered to first order (i.e. inter-facial effects are ignored). The properties of water are dependent on temperature, and this dependence is taken into account when computing w . The term 0 has the following dependence on temperature T (Hasted, 1973) 0 ðTÞ ¼ 87:74 0:40008T þ 9:398 104 T 2 1:41 106 T 3 ,
(13)
where the temperature is in Celsius. The dependence of the temperature on 1 and on q can be found in Table 2.2 of Hasted (1973). Besides, the influence of concentration on the permittivity of water is relatively small compared with the influence on the conductivity (Scho¨n, 2004), and only the effect on 0 is considered at this point. The correction term is (Olhoeft, 1981) 0 ¼ 0 ðTÞ 13s þ 1:065s2 0:03006s3 ,
porosity through a dimensionless number fL2 , (15) Tk which can be assumed constant for a given material (Johnson et al., 1986; Charlaix et al., 1988). Charlaix et al. (1988) obtained experimentally values of m 8 for fused glass beads, and m 10 for crushed glass. A value of m ¼ 8 can be obtained theoretically for a cylindrical tube (Pride, 1994). We have added the possibility to save the hydraulic permeability field obtained from Eq. (15) and the simulated porosity field, for user defined values of T and m. The permeability can also be expressed in terms of specific surface, porosity and tortuosity through the well-known Kozeny–Carman (K–C) relationship. The relation is in this case (Scho¨n, 2004) m¼
k¼
(16)
where w is the so-called Kozeny–Carman shape factor. w is equal to 2 for a circular cross-section (see Scho¨n, 2004 for values of w for other simple shapes). We have compared the permeabilities obtained with Eqs. (15) and (16) for a clean sand model (porosity field of Fig. 2a). The models are shown in Fig. 4. There is roughly a factor of 2.5 between the two, with lower values for the K–C model. This was observed by Raffensperger and Ferrell (1991), who applied the K–C model to experimental data and found predicted permeabilities 2–3 orders of magnitude lower than measured Salinity [ppm] 0
3900 3400
5
2900
(14)
where the concentration s is in mol/l.
f3 , wð1 fÞ2 S2 T2
10 Z [m]
¼ 1 þ
1087
2400 15 1900
2.4. Hydraulic permeability model 20
In his development, Pride (1994) used volume averaging to define three parameters related to the pore-space geometry: the ‘‘L parameter’’, the tortuosity T and the permeability k. His definitions are consistent with the previous work of Johnson et al. (1986, 1987). All three parameters are linked to
1400
25
900 0
5
10
15
20
25
X [m] Fig. 6. Example of user defined salinity model.
ARTICLE IN PRESS B. Giroux, M. Chouteau / Computers & Geosciences 34 (2008) 1080–1092
1088
σ [S/m]
σ [S/m]
0
0
0.020
0.07
0.019 5
0.06
5
0.018
0.05
0.017 Z [m]
0.016 0.015
15
10 Z [m]
10
0.04 15 0.03
0.014 0.013
20
20
0.02
0.012 25
25
0.011 0
5
10
15
20
25
0.01 0
5
10
X [m]
15
20
25
X [m]
Re(ε)
Re(ε) 11.0
0
11.0
0
10.5
10.5 5
5
10.0 9.5 9.0
15
9.0 15
8.5 8.0
20
9.5
10 Z [m]
Z [m]
10
10.0
8.5 8.0
20
7.5
7.5 25
25
7.0 0
5
10
15
20
25
7.0 0
5
10
X [m]
15
20
25
X [m] Im(ε)
Im(ε) 0.90
0
0.90
5
0.75
5
0.75
10
0.60
10
0.60
15
0.45
15
0.45
20
0.30
20
0.30
0.15
25
Z [m]
Z [m]
0
25 0
5
10
15 X [m]
20
25
0.15 0
5
10
15
20
25
X [m]
Fig. 7. Influence of salinity on electrical conductivity and relative dielectric permittivity: (a), (c), (e) Constant salinity at 900 ppm and (b), (d), (f) Salinity model of Fig. 6.
ARTICLE IN PRESS B. Giroux, M. Chouteau / Computers & Geosciences 34 (2008) 1080–1092
3. Examples 3.1. Influence of clay content We have presented in Section 2.1.1 how clay could be incorporated in the model, and how the clay fraction could be translated in terms of specific surface in Section 2.2. By default, specific surface values of 2:65 105 and 5 108 are used for Ss and Sc , respectively (Santamarina et al., 2002). It is, however, possible to model the conductivity of a clean sand aquifer by giving a value of Sc equal to that of Ss . In this section, the effect of clay is illustrated by comparing the conductivity and permittivity models obtained in clean sand and with the above default values of specific surface. The porosity model of Fig. 2a and clay fraction model derived from Eq. (2) (Fig. 3) are used for this example. The Cole–Cole parameters for the relative permittivity of clay are as following: 0 ¼ 30, 1 ¼ 20, t ¼ 108 and q ¼ 0:9. The electrolyte is composed of water and Naþ and Cl ions at 900 ppm concentration, at a temperature of 15 1C. The porosity field is fully saturated. The electrical conductivity fields obtained are shown in Figs. 5a and b. The presence of clay has for effect to increase the conductivity in low porosity zones, where the clay content is higher. The overall trend is then contrary to what is seen for clean sands. The relative dielectric permittivity models are shown in Figs. 5c–f. The high values of 0 and 1 for clay has for effect to increase the overall permittivity values, as seen in the figures. There is also a stronger increase of the imaginary part owing to the contribution of the increased conductivity. 3.2. Influence of salinity Let us consider next the effect of variation of salinity. Recall that salinity influences the electric conductivity (Eq. (4)) through the bulk ionic concentration N, as well as the dielectric permittivity according to Eq. (14). By default in our program the salinity is constant for all cells. However, we have included the
possibility to define a salinity model that is independent of the porosity field. This feature can be used to simulate approximately models at different times during a saline tracer infiltration test. The salinity model is kriged from points provided by the user. Fig. 6 depicts a hypothetical salinity model that is used to illustrate the influence of salinity on the conductivity and permittivity models. Consider the porosity field shown in Fig. 2a, on which the application of a constant salinity model (900 ppm) is compared to the application of the model of Fig. 6. In both cases, the porosity field is fully saturated with an electrolyte composed of water and Naþ and Cl ions, at 15 1C. The effect of clay is not considered here. The conductivity and permittivity fields obtained are shown in Fig. 7. The effect of increased salinity on conductivity and on the imaginary part of permittivity is very strong as expected. A slight increase of the real part of permittivity is also observed in Fig. 7d.
3.3. Influence of temperature and saturation As a third example, the effects of temperature and saturation are illustrated. The clayey sand, low salinity model defined in paragraph 3.1 is compared to a model with exactly the same parameters except 0
1
Vadose zone
2 Z [m]
values. They attributed the difference to the presence of dead-end pores and the importance of flow through large conducting channels. Chapuis and Aubertin (2003) discuss in length the influence of the specific surface on permeability estimation through the K–C relation.
1089
3
4
Capillary zone
5 0
0.2
0.4 0.6 Saturation
0.8
1
Fig. 8. Example of user defined saturation profile model. Defining parameters are minimum and maximum values of saturation, capillary fringe thickness and water table depth.
ARTICLE IN PRESS B. Giroux, M. Chouteau / Computers & Geosciences 34 (2008) 1080–1092
1090
σ25°C - σ15°C [S/m]
σ [S/m] 0
0
0.038
0.007 0.003
5
5 0.030
-0.001 10 0.022
Z [m]
Z [m]
10
-0.005 15
15
-0.009 0.014
20
25
20
5
10
15 X [m]
20
-0.017
25
0.006 0
-0.013
25
0
5
10
15 X [m]
20
25
Re(ε)25°C - Re(ε)15°C
Re(ε) 0
17.5
0
5
15.5
5
10
13.5
10
15
11.5
0 -2
Z [m]
Z [m]
-4 -6 15 -8 20
25 0
5
10
15 X [m]
20
9.5
20
7.5
25
25
-10 -12 0
5
Im(ε)
10
15 X [m]
20
25
Im(ε)25°C - Im(ε)15°C
0
1.5
0
7.2
1.0 5
0.5
5
5.9
0.0 -0.5
10 4.6
-1.0
Z [m]
Z [m]
10
-1.5 15
15 3.3
20
-2.0 -2.5
20
-3.0 -3.5
25
2.0 0
5
10
15 X [m]
20
25
25
-4.0 0
5
10
15 X [m]
20
25
Fig. 9. Influence of temperature on electrical conductivity and relative dielectric permittivity. Model is defined in Section 3.1 (clayey sand, low salinity), with an increase of temperature from 15 to 25 1C and with an added saturation profile (see text for details). (a), (c), (e) 25 1C and (b), (d), (f) difference.
ARTICLE IN PRESS B. Giroux, M. Chouteau / Computers & Geosciences 34 (2008) 1080–1092
for temperature which is increased from 15 to 25 1C. A simple saturation profile is also applied to the porosity field of the high temperature model. The shape of this saturation profile is defined by three parameters (residual saturation, thickness of capillary fringe and depth of water table) and follows the shape of an error function in the capillary zone, as shown in Fig. 8. For illustrative purpose, the thickness of the capillary fringe was set to exaggerated value of 3 m in this example. The water table is in this case at 5 m, and the values of saturation in the vadose zone and below the water table are 0.3 and 1, respectively. The model of van Genuchten (1980) is also implemented in HGmod. In this model, the saturation profile follows the following function: SðhÞ ¼ ½1 þ ðahÞn m ,
1091
Acknowledgments The authors acknowledge the financial support of NSERC. The first author is grateful to E. Gloaguen for fruitful discussions on geostatistical simulation. Appendix A. Computation of C em and C os For a fluid with L ionic species, the excess conductance associated with the electro-migration of double-layer ions is (Pride, 1994) L X ezl z 2 ðezl Þ bl Nl exp C em ¼ 2d 1 , 2kT l¼1 (A.1) where
(17)
where h is the pressure head and a, m and n are empirical parameters. The results for 15 1C are shown above, in Figs. 5b, d and f (for fully saturated soils). Fig. 9 shows the conductivity and permittivity for 25 1C, as well as the differences between the 25 and 15 1C results. As shown in Fig. 9, the effect of increased temperature is to slightly increase the conductivity and imaginary part of permittivity in the saturated zone. The real part of permittivity is virtually unchanged for this 10 1C difference in temperature. In the unsaturated zone, the reduced water content has a marked effect on conductivity and permittivity, leading to a decrease of both, as expected.
4. Conclusion We have presented a computer program that permits to build a 2D or 3D grid model of dielectric permittivity, electrical resistivity and hydraulic permeability, all derived from a stochastic realization of a porosity field. The program is written in C++ and the source code is available at http:geo.polmtl.ca/giroux/HGmod. The program is relatively simple and modular, and can easily accommodate extensions. Anticipated future additions are the inclusion of water flow modeling capabilities to generate snapshots of the properties at various times that can simulate an infiltration experiments, for example. We also investigate to include dispersive complex conductivity so as to model the subsurface response to induced polarization.
d is the Debye length, which is a measure of the thickness of the diffuse double-layer and which is defined as L X 1 ðezl Þ2 Nl , ¼ k kT d2 l¼1 0 f
(A.2)
with 0 the dielectric permittivity of vacuum, kf the dielectric constant of the fluid, k Boltzmann’s constant ð1:38 1023 J=KÞ, and T the temperature in Kelvin; z is the zeta potential (in Volt), for which the expression 0:008 þ 0:026 log10 s can be used (Carcione et al., 2003), in which s is salt concentration in mol/l.
For the same solution, the conductance due to electrically induced streaming (convection) of the excess double-layer ions (the electro-osmotic conductance) is 1 ð0 kf zÞ2 P 2i3=2 d C os ¼ 1 , (A.3) Pd 2dZf where
P is a simplification term equal to 8kTd 2 X ezl z P¼ LNl exp 1 , (A.4) 2kT 0 kf z2 l¼1
d is the viscous skin depth, equal to sffiffiffiffiffiffiffiffiffi Zf d¼ , orf
(A.5)
ARTICLE IN PRESS 1092
B. Giroux, M. Chouteau / Computers & Geosciences 34 (2008) 1080–1092
with o the angular frequency and rf the fluid density. References Archie, G.E, 1942. The electrical resistivity log as an aid in determining some reservoir characteristics. Transactions of American Institute of Mining Metallurgical Engineers 146, 54–62. Carcione, J.M., Seriani, G., Gei, D., 2003. Acoustic and electromagnetic properties of soils saturated with salt water and NAPL. Journal of Applied Geophysics 52, 177–191. Chang, P.-Y., Alumbaugh, D., Brainard, J., Hall, L., 2006. Cross-borehole ground-penetrating radar for monitoring and imaging solute transport within the vadose zone. Water Resources Research 42, W10413. Chapuis, R.P., Aubertin, M., 2003. On the use of the Kozeny–Carman equation to predict the hydraulic conductivity of soils. Canadian Geotechnical Journal 40 (3), 616–628. Charlaix, E., Kushnick, A.P., Stokes, J.P., 1988. Experimental study of dynamic permeability in porous media. Physical Review Letters 61 (14), 1595–1598. Chile`s, J.-P., Delfiner, P., 1999. Geostatistics: Modeling Spatial Uncertainty. Wiley, New York, 695 pp. Day-Lewis, F.D., Singha, K., Binley, A.M., 2005. Applying petrophysical models to radar travel time and electrical resistivity tomograms: resolution-dependent limitations. Journal of Geophysical Research 110, B08206. de Lima, O.A.L., Sharma, M.M., 1990. A grain conductivity approach to shaly sandstones. Geophysics 55 (10), 1347–1356. Druchinin, S.V., 2000. Models for calculation of dielectric constant of moist sandy–clayey soils in wavelengths from centimeters to tens of meters. In: Noon, D., Stickley, G., Longstaff, D. (Eds.), Proceedings of the Eighth International Conference on GroundPenetrating Radar, Gold Coast, Australia. Freeze, R.A., 1975. A stochastic-conceptual analysis of onedimensional groundwater flow in nonuniform homogeneous media. Water Resources Research 11 (5), 725–741. Furman, A., Ferre´a, T.P.A., Warrick, A.W., 2004. Optimization of ERT surveys for monitoring transient hydrological events using perturbation sensitivity and genetic algorithms. Vadose Zone Journal 3, 1230–1239. Hasted, J.B., 1973. Aqueous Dielectrics. Chapman & Hall, London, 302 pp. Hubbard, S., Rubin, Y., Majer, E., 1999. Spatial correlation structure estimation using geophysical data. Water Resources Research 35 (6), 1809–1825. Johnson, D.L., Sen, P.N., 1988. Dependence of the conductivity of a porous medium on electrolyte conductivity. Physical Review B 37 (7), 3502–3510. Johnson, D.L., Koplik, J., Schwartz, L.M., 1986. New pore-size parameter characterizing transport in porous media. Physical Review Letters 57 (20), 2564–2567. Johnson, D.L., Koplik, J., Dashen, R., 1987. Theory of dynamic permeability and tortuosity in fluid-saturated porous media. Journal of Fluid Mechanics 176, 379–402.
Le Ravalec, M., Noetinger, B., Hu, L.Y., 2000. The FFT moving average (FFT-MA) generator: an efficient numerical method for generating and conditioning Gaussian simulations. Mathematical Geology 32 (6), 701–723. Marion, D., Nur, A., Yin, H., Han, D., 1992. Compressional velocity and porosity in sand–clay mixtures. Geophysics 57 (4), 554–563. Olhoeft, G.R., 1981. Electrical properties of rocks. In: Touloukian, Y.S., Judd, Y.S., Roy, R.F. (Eds.), Physical Properties of Rocks and Minerals. McGraw-Hill, New York, pp. 257–330. Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., 1992. Numerical Recipes in C: The Art of Scientific Computing, second ed. Cambridge University Press, Cambridge, 994 pp. Pride, S., 1994. Governing equations for the coupled electromagnetics and acoustic of porous media. Physical Review B 50 (21), 15678–15696. Pride, S.R., Morgan, F.D., Gangi, A.F., 1993. Drag forces of porous-medium acoustics. Physical Review B 47 (9), 4964–4978. Raffensperger, J.P., Ferrell Jr., R.E., 1991. An empirical model of intrinsic permeability in reactive clay-bearing sands. Water Resources Research 27 (11), 2835–2844. Revil, A., 1999. Ionic diffusivity, electrical conductivity, membrane and thermoelectric potentials in colloids and granular porous media: a unified model. Journal of Colloid and Interface Science 212 (2), 503–522. Santamarina, J., Klein, K., Wang, Y., Prencke, E., 2002. Specific surface: determination and relevance. Canadian Geotechnical Journal 39, 233–241. Scho¨n, J.H., 2004. Physical properties of rocks: fundamentals and principles of petrophysics, Handbook of Geophysical Exploration—Seismic Exploration, vol. 18, first ed. Elsevier, Amsterdam, 583 pp. Sihvola, A., 2000. Mixing rules with complex dielectric coefficients. Subsurface Sensing Technologies and Applications 1 (4), 393–415. Sun, H., Koch, M., 1998. Fractal generation of surface area of porous media. Stochastic Hydrology and Hydraulics 12, 83–96. Tabbagh, A., Panissod, C., Gue´rin, R., Cosenza, P., 2002. Numerical modeling of the role of water and clay in soils’ and rocks’ bulk electrical conductivity. Journal of Geophysical Research 107 (B11), 2318. Tronicke, J., Holliger, K., 2005. Quantitative integration of hydrogeophysical data: conditional geostatistical simulation for characterizing heterogeneous alluvial aquifers. Geophysics 70 (3), H1–H10. van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal 44, 892–898. Waxman, M.H., Smits, L.J.M., 1968. Electrical conductivities in oil-bearing shaly sands. Society of Petroleum Engineering Journal 8, 107–122. Weast, R., Astle, M.J. (Eds.) 1980. Handbook of Chemistry and Physics. 61st ed., CRC Press, Boca Raton, Florida, 235 pp.
