- Email: [email protected]
Pore-scale modeling of electrical conductivity in unsaturated sandstones
235
P o r e - s c a l e m o d e l i n g of e l e c t r i c a l c o n d u c t i v i t y in u n s a t u r a t e d
sandstones
G. Cassiani ~, E. Dalla b, A. Brovelli b, and D. Pitea b* ~Dipartimento di Scienze Geologiche e Geotecnologie, Univ. di Milano-Bicocca Piazza della Scienza 4, 20126 Milano, Italy bDipartimento di Scienze dell'Ambiente e del Territorio, Univ. di Milano-Bicocca P.za della Scienza 1, 20126 Milano Non-invasive techniques are becoming common tools for site characterization. These techniques include electrical and electromagnetic methods, and are generally referred to as hydro-geophysical techniques. This approach relies upon the correct translation of geophysical properties into hydrological ones, and particularly moisture content (in the vadose zone) and solute concentration (in both the vadose and saturated zones). In this note, we present a pore-scale approach for the investigation of the relationships that link electrical properties of porous media with water content and the solid interfacial areas responsible for surface conductivity. In this approach we create a digital porous medium on the basis of experimental grain size distribution data and apply a previously developed drainage simulator based on pore-morphology analysis. Then we simulate electrical flow within the pore-scale distributions of the phases and derive the bulk resistivity of the porous medium. Simulation results are compared with resistivity data from a UK PermoTriassic Sandstone (Sherwood S.); resistivity of the solid phase, or alternatively surface conductivity, has been derived from calibration against measured bulk resistivity. The agreement between measured and simulated bulk conductivities, at different saturation values, is very good. 1. I N T R O D U C T I O N
The use of geophysical methods for site characterization has a long tradition in the hydrological community [1]. However, only the recent development of powerful field equipment, processing, and inversion methods has brought the techniques up to the requirement for a detailed, non-invasive imaging of the subsurface in its static structure (geology) and dynamic processes (hydrology). Among such new techniques, electrical resistivity tomography (ERT), both in surface and cross-hole configuration, has gained respect in the hydrological community as a tool to investigate both the unsaturated and the saturated *We acknowledge support from A. Binley (Lancaster University), who helped define the case study. Funding from the UK Natural Environment Research Council Grant NER/A/S/2001/00411 was partly devoted to the work presented herein. This work was also supported by the Consorzio Interuniversitario CINECA (Bo- Italy). We thank N. Fusi and G. Crosta for their help on the experimental measurement of capillary pressure-saturation curve.
236 zones. For example, above the water table, ERT is sensible to the water content distribution, and can be used to study the dynamics of infiltration [2-7]. The use of non invasive techniques is attractive, because they allow for a spatially extensive characterization of the subsurface, generally at the scale of interest (meters to tens of meters), with limited disturbance to the natural processes. However, such new applications are possible only if a suitable constitutive model can be used to translate the geophysical property (for ERT, DC resistivity) into a hydrological property (e.g., moisture content). Traditionally, the dependence of bulk resistivity upon porous medium characteristics, porosity, and water saturation has been expressed by the empirical Archie's law [8]. Archie's law expresses the bulk conductivity, orb, of a fully saturated sample (water saturation s~ - 1) as" O"w
Ob(S~- 1 ) - F
(1)
where cr~ is the conductivity of the aqueous phase and F (dimensionless) is the electrical formation factor. F is generally expressed as F - ,1, where ~ is the porosity and m is named the cementation exponent. The cementation exponent is known to be a function of the connectivity of the conductive phase within the porous medium. Archie's law can be extended to partially saturated media (Sw < 1) by defining the resistivity index RI as: Orb ir~I--
O-b(S w - -
1) - s ; ~
(2)
where n is an empirical parameter usually referred to as the Archie saturation exponent (related to the connectivity of the porous medium), with values close to 2 in most cases. The underlying assumption of Archie's law is that only one conducting fluid phase exists distributed within a non-conducting matrix. This approximation is not acceptable for the study of shaly sands, where a significant fraction of clay minerals is present. In this case the conduction pathway due to the electrical double layer formed at the solid phase-aqueous phase interface is significant [9]. Different models have been proposed in the literature, such as the classical Waxman and Smith model [10], the dual water-model [11], variations of these models [4,12], or models based upon differential theory [13]. The objectives of this paper are: 1. to analyze the electrical conductance in a porous medium starting from basic principles and simple information (porosity, grain size distribution), by building a model that represents the medium at the pore scale; 2. to analyze the effect and relative importance of solid conductivity, modeled either as solid bulk conductivity, or as surface conductivity; 3. to use the model to predict the change of bulk conductivity as a function of water saturation, in order to compare the predicted behavior with classical models. The tool of choice in order to achieve the above objectives is a pore-scale model based on a random packing of spheres [14] and on the fine discretization of the pore space, the solid matrix, and the solid-fluid interfaces. This type of pore-scale model has been recently recognized in the hydrogeology literature as a useful tool for bridging the gap
237 between microscale phenomena and macroscale description of the system, with its bulk properties. To date, only few studies [15-19] investigated electrical properties of a porous medium with a pore-scale approach. These studies do not include any comparison with experimental data. In this paper, simulation results are compared with experimental measurements from weakly consolidated Permo-Triassic sandstones from the UK, with a particular attention toward calibrating the electrical characteristics of the solid phase. 2. M E T H O D O L O G Y
We adopted a five-step approach to model electrical conductivity as a function of water saturation in the porous medium: 1. We created a synthetic, non-overlapping and gravitationally stable packing of spheres using a random packing code [14], honoring prescribed porosity and grain size distribution. This virtual porous medium is then discretized into elementary cubic volumes ("voxels") using a code that preserves the original porosity. .
We simulated primary drainage by applying a simulator [20] that performs a morphological analysis of the digital pore space. This approach assumes the bottom of the digital porous medium to be connected to a non wetting phase (NWP) reservoir and the top to a wetting phase (WP) reservoir. At the beginning of the simulation we assume the same pressure exists in both reservoirs, PNWP PwP, and the porous medium saturated with the wetting phase. As the pressure in the non wetting phase reservoir is increased, and thus the capillary pressure, p~ = P N W P - PWF, the non wetting phase invades the domain. The new pore-scale distribution of the phases is computed using Laplace equation for a spherical interface and zero contact angle. The main assumptions of the code are that (a) the porous medium is totally water wet and (b) no irreducible wetting phase is present, s~ - 0. The simulated primary drainage curve was compared with experimental results on samples of chosen material.
3. We developed a steady-state simulator that models bulk electrical current through the porous medium. This model is based upon a cell-centered finite difference discretization of the 3D flow equation (V is the electrical potential):
v.
(o-vv) - o
(3)
where the volume electrical conductivity, a, expressed in Siemens per meter, at the voxel interfaces is computed as the harmonic mean of the conductivity of two adjacent voxels. We assume phase boundaries to be voxel boundaries and compute the potential at the voxel centers. Dirichlet boundary conditions are applied to two opposite faces of the porous medium, and no-flow boundaries are prescribed on the remaining four faces. Solution of the system of equations is obtained by using a projection method with LU preconditioner.
238
o
We developed a steady-state simulator that models electrical current flowing along the solid-liquid interfaces within the porous medium. The surface representing the interfacial area between solid phase and aqueous phase is represented as the triangular mesh generated by the standard marching-cube algorithm [21] applied to the porous medium already discretized in voxels. The method utilized to compute interfacial areas in digital porous media has been fully validated [22]. The resulting structure is a triangular irregular 3D network, suitable for numerical simulation via a Galerkin Finite Element formulation of the equation V . (o~VV) - 0
(4)
where cri is the surface conductivity per unit length (expressed in Siemens), and the nabla operators are applied in 2D along the surfaces of two-dimensional triangular linear elements connected in the 3D space. Dirichlet boundary conditions are applied to the nodes lying on two opposite faces of the porous medium, and no flow boundary conditions are applied to the remaining four faces. 5. From the simulation results we compute the net flux of electrical-charge moving through the digital porous medium. The volumes contribution to the bulk conductivity is then simply derived, according to Ohm's law, as I
1
= AV s
(5)
where I is the current intensity, A V is the potential difference applied, 1 is the length of the domain, and S is the area of each surface where a first-type (Dirichlet) boundary condition is applied. The contribution to the bulk conductivity of the conductance along the surfaces is computed by using the same equation 5 where I is now the total electrical current flowing in or out of the triangulated surface network at the locations of the Dirichlet boundary conditions. 3. S I M U L A T I O N RESULTS" P O R O U S M E D I U M A N D P R I M A R Y D R A I N A G E
As a test material, we selected the UK Permo-Triassic Sandstones (Sherwood S.), which constitute the second largest aquifer system in the UK. This type of porous medium is representative of a fine to medium, poorly consolidated sandstone, with a small fraction of clay and Fe oxide material coating the silicate grains. Consequently, the Sherwood Sandstone is expected to exhibit a non-null solid phase/surface conductivity, suitable for analysis with the model described in Section 2. Samples of the Sherwood Sandstone from South Yorkshire [23] were analyzed for porosity and grain-size distribution. We assumed a log-normal distribution for the grain size distribution and, thus, simulated a gravitationally stable sphere packing that follows the logarithmic mean, standard deviation, and porosity. Simulated packing properties are listed in Table 1. The primary drainage curve was experimentally measured on three remolded Sherwood sandstone samples with a tempe pressure cell apparatus. The non-wetting phase was air
239 Table 1 Properties of the simulated sphere-packing. Par amet er Value U nit grain radius Rg 0.154 + 0.051 mm domain length L 7.974 mm porosity 9 0.39 number of spheres 14501
and the wetting phase was distilled water. The surface tension between air and water, "7 has been assumed equal to 72 x 10 -a N/m. The porous medium described in Table 1 was discretized with increasing resolution of 1283, 2563 and 3843 voxels to simulate the primary drainage curve. Figure 1 compares experimental results with simulated curves at three discretization levels. Results are affected by resolution only at very low water saturation, while the overall shape of the drainage curve is not substantially affected by the discretization level. Analogous discretization effects on simulated p ~ - s ~ curves were observed and discussed in the literature [20]. Our simulations reproduce correctly the features of the experimental data. Particularly, we obtain an excellent agreement for the flat portion of the pC - s ~ curve and for simulated and experimental air entry pressure value, which are related to the uniformity of the medium and to the mean pore throat diameter, respectively. This implies that the simulated packing we used is well representative of the real pore structure.
9
o O <) 9 9 V__
~-250m-: 3:~ / ] E
~aooJ: ~....
9
Simulated data, 1283 voxels Simulated data, 2563 voxels Simulated data, 3843 voxels Experimental data, 1st set Experimental data, 2 nd set
Experi___ment_aldata__,3:dset_j__
lS~ a,IO%
~
.,,
/
~+'~' ~ .o,~+o~,.,~,"
%
o'.a
o'.4
o'.+
Water saturation
o'.a
-'r
Figure 1. Discretization effects on primary drainage curve and comparison with three sets of experimental data.
240 4. S I M U L A T I O N
RESULTS" BULK ELECTRICAL
CONDUCTIVITY
We computed the bulk electrical conductivity of the porous medium for each equilibrium air-water distribution modeled with the morphological simulator. The three-dimensional simulator discussed in Section 2 was set to run on an IBM SP RS/6000 Power4. 4.1. C a l i b r a t i o n of m a t r i x e l e c t r i c a l c o n d u c t i v i t y Two approaches are possible for the description of conductance related to the clay/oxide fraction of the porous medium" 1. surface conductivity is represented by an equivalent conductivity of the whole solid matrix volume, a~; 2. surface conductivity is represented as a proper interfacial property along the solid/water interfaces, cq. We investigated both approaches. In either cases, the unknown parameter needs to be determined by matching the experimental value of saturated bulk conductivity of the sandstone [23], equal to a b - 1.56 x 10 .2 S/re. 4.1.1. solid e q u i v a l e n t c o n d u c t i v i t y a p p r o a c h This calibration was performed on a completely saturated porous medium, discretized with 3843 voxels. We assume aw - 6.8 x 10 -2 S//m, as from experimental measurements by [23] for natural groundwater at the site. Figure 2a shows the results. The resulting bulk solid conductivity is a~ - 1.43 x 10 -3 S/re. As expected, this is a nonzero value, even though significantly lower than water conductivity. In this approach, ~ represents the additional surface conductivity due to mobile ions in the electrical double layer of the clay.
0.019
0.016-
0.0158 ~
0.018
,,,,~' ,,
0.0156 . . . . . . . . . . . . . . . . . . . . . . . .
,,,
0.017
,,'0'"
,
,.-- 0.0154 E
'~'~.016 _
,,,"
,,"'"
"-".,,0.0152 "5 0.015
,,'"
0.015 ~ 0.0148-
0.014
,,,
,,"
,,~'
,,,'" O"
0.0146-
a)
o.o~3o
0'.5
1
,I', 1.5 ~solid (S/m)
2.5
3
x 0 b)
5
6
7
o,s
x 04 '
Figure 2. Calibration of solid phase conductivity for a fully saturated medium with the solid equivalent conductivity approach (a) and the interface conductivity approach (b).
241 Different models have been proposed to take into account the electrical double layer conduction. The effective bulk conductivity of a saturated medium can be defined [10] as
1
~v - T ( B Q v + ~ )
(6)
where Q~ is the volume concentration of clay exchange cations (meq/mL) and B the equivalent conductance of the clay counterions (cm2S/meq), which is related to their mobility. Some works [12] define B Q ~ - (Tclay and compute it assuming conduction pathways through the solid phase. According to this model, our cr~ - cr~lay. The cr~ value we compute is lower than those presented in both papers, indicating a really low surface conductance. For simulations of electrical charge flow under unsaturated conditions, we assumed cr~lay = BQ~ ~ [10]. This states that when water saturation is less than unit, the exchange cations associated with the clay become more concentrated in the remaining pore water.
4.1.2. interface conductivity approach We made also a first attempt at calibrating the surface conductivity according to the second approach, on a small porous medium derived as a 15a voxels subdomain of the digital porous medium discretized with 128 voxels per side. Figure 3 shows the triangular mesh obtained with the marching cubes code that models the interface between solid and aqueous phase.
Figure 3. Triangular mesh used to compute the surface conductivity.
242
Calibration results are shown in Figure 2b. We calculated ai - 9.0 x 10 . 4 S. This value is of the same order of magnitude of the one calculated for a volume equivalent conductivity of the solid phase. Due to the preliminary nature of the results a more specific comparison between the two values may not be meaningful. Particularly investigation of domain size effects on simulated ai deserves further studies. Thus, at this stage, comparison against experimental data has been performed using approach 1. 4.2. F i c t i t i o u s e l e c t r i c a l c o n d u c t i v i t y of g a s e o u s p h a s e An effective value for electrical conductivity of the gaseous phase (air) has to be determined for numerical reasons, since a value equal to zero would possibly result in a singular matrix of the finite difference system of equations. We adopted a~ - 6.8 • 10 -2 S / m and a~ - 1.43 • 10 -3 S/m and calibrated the value of gaseous phase electrical conductivity on a porous medium discretized with 3843 voxels with s~ - 0.13 (a high air content). The gas conductivity, physically zero, must be chosen as large as possible in order not to generate a matrix too ill-conditioned, but must be small enough so that any smaller gas conductivity would effectively yield the same medium bulk conductivity (asymptotic behavior) (see Figure 4). The asymptote is reached for O"a nearly equal to 1.0 • 10 -7 S/m, typical of an insulator. Note that the computational time increases with decreasing gaseous phase conductivity~ as the system matrix condition number increases, and forces to the iterative solver to converge more slowly.
x 10 -4 ,
i
9.54 9.52
....-
9.5
.-
9.48 E (/)
9.46
,,,,,,,'~
9.44 ,3
9.42 9.4 9.38 9.36 10 -7
lO-(s/m )
10 -5
(~air
Figure 4. Calibration of gaseous phase conductivity.
4.3. C o m p a r i s o n w i t h e x p e r i m e n t a l d a t a The results of electrical field simulations on the 3843 voxels porous medium are shown in Figure 5 as a function of water saturation~ together with experimental measurement
243 results [23]. Comparison between simulated and experimental bulk conductivities show an excellent agreement at high and intermediate s ~ values.
450 400
I
- O " Simulated data 9 Experimental data
350 300 E
~250 -5 .o (:2.
200 150 100 1
0.4 0.6 Water saturation
i
0.8
1
Figure 5. Simulated bulk conductivities and experimental data from [23] as a function of water saturation.
5. C O N C L U S I O N S
A pore-scale method based on digital representation of a porous medium was developed and applied to predict the bulk electrical conductivity of unsaturated porous media. The model was used to simulate the DC resistivity of semi-consolidated Triassic sandstone from the UK (Sherwood sandstone). Comparison with experimental data show a very good agreement with model results. The results presented in this paper show the accuracy and reliability of the proposed method for pore-scale modeling and computation of bulk conductivities of a porous medium. This method can be a powerful tool in the investigation of constitutive electrical law from a pore scale point of view. Further work will focus on testing the surface conductivity model and extend the model to account for induced polarization phenomena. REFERENCES
1. O. Mazac, W. E. Kelly, I. Landa Surface geoelectrics for groundwater pollution and protection studies Journal of Hydrology, vol.93, no.3-4, pages.277-294, 1987. 2. R.K. Frohlich and C. D. Parke. The electrical resistivity of the vadose zone - field survey. Ground Water, 27(4):524-530, 1989. 3. W. D. Daily, A. L. Ramirez, D. J. LaBrecque, and J. Nitao. Electrical resistivity tomography of vadose water movement. Water Resour. Res., 28:1429-1442, 1992.
244 4.
R. J. Kalinski and W. E. Kelly. Estimating water content of soils from electrical resistivity. Geotechnical testing journal, pages 323-329, 1993. 5. A. Ramirez, W. Daily, A. Binley, D. LaBrecque, and D. Roelant. Detection of leaks in underground storage tanks using electrical resistence methods. J. Env. and Eng. Geophysics, 1(3)'189-203, 1996. 6. A. Binley, G. Cassiani, R. Middleton, and P. Winship. Vadose zone model parameterisation using cross-borehole radar and resistivity imaging. J. Hydrol, 267(3-4)" 147-159, 2002. 7. D. Michot, Y. Benderitter, A. Dorigny, B. Nicoullaud, D. King, and A. Tabbagh. Spatial and temporal monitoring of soil water content with an irrigated corn crop cover using surface electrical resistivity tomography. Water Resources Research, 39
(5), 2003. 8. 9. 10. 11.
12. 13.
14. 15.
16.
17. 18. 19. 20. 21. 22.
G. E. Archie. The electrical resistivity log as an aid in determining some reservoir characteristics. Trans. Am. Inst. Min. Metall. Pet. Eng., 146:54-62, 1942. A.E. Bussian. Electrical conductance in a porous medium. Geophysics, 48:1258-1268, 1983. M. H. Waxman and L. J. M. Smits. Electrical conductivities in oil-bearing shaly sands. Soc. Pet. Eng. J., 243:107-122, 1968. C. Clavier, G. Coates, and J. Dumanoir. The theoretical and experimental bases fro the dual-water model for the interpretation of shaly sands. Soc. Pet. Eng. J., 24(2)" 153-169, 1984. P. W. J. Glover, M. J. Hole, and J. Pous. A modified archie's law for two conducting phases. Earth and planetary sciences letters, 180:369-383, 2000. A. Revil, L. M. Cathles III, and S. Losh. Electrical conductivity in shaly sands with geophysical applications. Journal of geophysical research, 103(B10)'23,925-23,936, 1998. A. Yang, C. T. Miller, and L. D. Turcoliver. Simulation of correlated and uncorrelated packing of random size spheres. Physical Review E, 53(2)'1516-1524, 1996. A. L. Endres and R. Knight. The effects of pore-scale fluid distribution on the physical properties of partially saturated tight sandstones. J. Appl. Phys., 69(2)'1091-1098, 1991. R. J. Suman and R. J. Knight. Effects of pore structure and wettability on electrical resistivity of partially saturated rocks - a network study. Geophysics, 62(4)'1151-1162, 1997. D. Zhou, S. Arbabi, and E. H. Stenby. A percolation study of wettability effect on the electrical properties of reservoir rocks. Transport in porous media, 29:85-98, 1997. C. D. Tsakiroglou and M. Fleury. Pore network analysis of resistivity index for waterwet porous media. Transport in porous media, 35:89-128, 1999. H. N. Man and X. D. Jing. Pore network modelling of electrical resistivity and capillary pressure characteristics. Transport in Porous Media, 41:263-286, 2000. M. Hilpert and C. T. Miller. Pore-morphology-based simulation of drainage in totally wetting porous media. Advances in Water Resources, 24:243-255, 2001. Lorensen, W. E. and Cline, H. E. Marching cubes: A high resolution 3D surface construction algorithm. Computer Graphics, 21:163-169, 1987. E. Dalla, M. Hilpert, and C. T. Miller. Computation of the interracial area for multi-
245 fluid porous medium systems. Journal of Contaminant Hydrology, 56(1-2)'25-48, 2002. 23. A. Binley, P. Winship, R. Middleton, M. Pokar, and J. West. Observation of seasonal dynamics in the vadose zone using borehole radar and resistivity. In Proceedings
of the symposium on the application of geophysics to environmental and engineering problems, SAGEEP, march 4-7 2001.
P o r e - s c a l e m o d e l i n g of e l e c t r i c a l c o n d u c t i v i t y in u n s a t u r a t e d
sandstones
G. Cassiani ~, E. Dalla b, A. Brovelli b, and D. Pitea b* ~Dipartimento di Scienze Geologiche e Geotecnologie, Univ. di Milano-Bicocca Piazza della Scienza 4, 20126 Milano, Italy bDipartimento di Scienze dell'Ambiente e del Territorio, Univ. di Milano-Bicocca P.za della Scienza 1, 20126 Milano Non-invasive techniques are becoming common tools for site characterization. These techniques include electrical and electromagnetic methods, and are generally referred to as hydro-geophysical techniques. This approach relies upon the correct translation of geophysical properties into hydrological ones, and particularly moisture content (in the vadose zone) and solute concentration (in both the vadose and saturated zones). In this note, we present a pore-scale approach for the investigation of the relationships that link electrical properties of porous media with water content and the solid interfacial areas responsible for surface conductivity. In this approach we create a digital porous medium on the basis of experimental grain size distribution data and apply a previously developed drainage simulator based on pore-morphology analysis. Then we simulate electrical flow within the pore-scale distributions of the phases and derive the bulk resistivity of the porous medium. Simulation results are compared with resistivity data from a UK PermoTriassic Sandstone (Sherwood S.); resistivity of the solid phase, or alternatively surface conductivity, has been derived from calibration against measured bulk resistivity. The agreement between measured and simulated bulk conductivities, at different saturation values, is very good. 1. I N T R O D U C T I O N
The use of geophysical methods for site characterization has a long tradition in the hydrological community [1]. However, only the recent development of powerful field equipment, processing, and inversion methods has brought the techniques up to the requirement for a detailed, non-invasive imaging of the subsurface in its static structure (geology) and dynamic processes (hydrology). Among such new techniques, electrical resistivity tomography (ERT), both in surface and cross-hole configuration, has gained respect in the hydrological community as a tool to investigate both the unsaturated and the saturated *We acknowledge support from A. Binley (Lancaster University), who helped define the case study. Funding from the UK Natural Environment Research Council Grant NER/A/S/2001/00411 was partly devoted to the work presented herein. This work was also supported by the Consorzio Interuniversitario CINECA (Bo- Italy). We thank N. Fusi and G. Crosta for their help on the experimental measurement of capillary pressure-saturation curve.
236 zones. For example, above the water table, ERT is sensible to the water content distribution, and can be used to study the dynamics of infiltration [2-7]. The use of non invasive techniques is attractive, because they allow for a spatially extensive characterization of the subsurface, generally at the scale of interest (meters to tens of meters), with limited disturbance to the natural processes. However, such new applications are possible only if a suitable constitutive model can be used to translate the geophysical property (for ERT, DC resistivity) into a hydrological property (e.g., moisture content). Traditionally, the dependence of bulk resistivity upon porous medium characteristics, porosity, and water saturation has been expressed by the empirical Archie's law [8]. Archie's law expresses the bulk conductivity, orb, of a fully saturated sample (water saturation s~ - 1) as" O"w
Ob(S~- 1 ) - F
(1)
where cr~ is the conductivity of the aqueous phase and F (dimensionless) is the electrical formation factor. F is generally expressed as F - ,1, where ~ is the porosity and m is named the cementation exponent. The cementation exponent is known to be a function of the connectivity of the conductive phase within the porous medium. Archie's law can be extended to partially saturated media (Sw < 1) by defining the resistivity index RI as: Orb ir~I--
O-b(S w - -
1) - s ; ~
(2)
where n is an empirical parameter usually referred to as the Archie saturation exponent (related to the connectivity of the porous medium), with values close to 2 in most cases. The underlying assumption of Archie's law is that only one conducting fluid phase exists distributed within a non-conducting matrix. This approximation is not acceptable for the study of shaly sands, where a significant fraction of clay minerals is present. In this case the conduction pathway due to the electrical double layer formed at the solid phase-aqueous phase interface is significant [9]. Different models have been proposed in the literature, such as the classical Waxman and Smith model [10], the dual water-model [11], variations of these models [4,12], or models based upon differential theory [13]. The objectives of this paper are: 1. to analyze the electrical conductance in a porous medium starting from basic principles and simple information (porosity, grain size distribution), by building a model that represents the medium at the pore scale; 2. to analyze the effect and relative importance of solid conductivity, modeled either as solid bulk conductivity, or as surface conductivity; 3. to use the model to predict the change of bulk conductivity as a function of water saturation, in order to compare the predicted behavior with classical models. The tool of choice in order to achieve the above objectives is a pore-scale model based on a random packing of spheres [14] and on the fine discretization of the pore space, the solid matrix, and the solid-fluid interfaces. This type of pore-scale model has been recently recognized in the hydrogeology literature as a useful tool for bridging the gap
237 between microscale phenomena and macroscale description of the system, with its bulk properties. To date, only few studies [15-19] investigated electrical properties of a porous medium with a pore-scale approach. These studies do not include any comparison with experimental data. In this paper, simulation results are compared with experimental measurements from weakly consolidated Permo-Triassic sandstones from the UK, with a particular attention toward calibrating the electrical characteristics of the solid phase. 2. M E T H O D O L O G Y
We adopted a five-step approach to model electrical conductivity as a function of water saturation in the porous medium: 1. We created a synthetic, non-overlapping and gravitationally stable packing of spheres using a random packing code [14], honoring prescribed porosity and grain size distribution. This virtual porous medium is then discretized into elementary cubic volumes ("voxels") using a code that preserves the original porosity. .
We simulated primary drainage by applying a simulator [20] that performs a morphological analysis of the digital pore space. This approach assumes the bottom of the digital porous medium to be connected to a non wetting phase (NWP) reservoir and the top to a wetting phase (WP) reservoir. At the beginning of the simulation we assume the same pressure exists in both reservoirs, PNWP PwP, and the porous medium saturated with the wetting phase. As the pressure in the non wetting phase reservoir is increased, and thus the capillary pressure, p~ = P N W P - PWF, the non wetting phase invades the domain. The new pore-scale distribution of the phases is computed using Laplace equation for a spherical interface and zero contact angle. The main assumptions of the code are that (a) the porous medium is totally water wet and (b) no irreducible wetting phase is present, s~ - 0. The simulated primary drainage curve was compared with experimental results on samples of chosen material.
3. We developed a steady-state simulator that models bulk electrical current through the porous medium. This model is based upon a cell-centered finite difference discretization of the 3D flow equation (V is the electrical potential):
v.
(o-vv) - o
(3)
where the volume electrical conductivity, a, expressed in Siemens per meter, at the voxel interfaces is computed as the harmonic mean of the conductivity of two adjacent voxels. We assume phase boundaries to be voxel boundaries and compute the potential at the voxel centers. Dirichlet boundary conditions are applied to two opposite faces of the porous medium, and no-flow boundaries are prescribed on the remaining four faces. Solution of the system of equations is obtained by using a projection method with LU preconditioner.
238
o
We developed a steady-state simulator that models electrical current flowing along the solid-liquid interfaces within the porous medium. The surface representing the interfacial area between solid phase and aqueous phase is represented as the triangular mesh generated by the standard marching-cube algorithm [21] applied to the porous medium already discretized in voxels. The method utilized to compute interfacial areas in digital porous media has been fully validated [22]. The resulting structure is a triangular irregular 3D network, suitable for numerical simulation via a Galerkin Finite Element formulation of the equation V . (o~VV) - 0
(4)
where cri is the surface conductivity per unit length (expressed in Siemens), and the nabla operators are applied in 2D along the surfaces of two-dimensional triangular linear elements connected in the 3D space. Dirichlet boundary conditions are applied to the nodes lying on two opposite faces of the porous medium, and no flow boundary conditions are applied to the remaining four faces. 5. From the simulation results we compute the net flux of electrical-charge moving through the digital porous medium. The volumes contribution to the bulk conductivity is then simply derived, according to Ohm's law, as I
1
= AV s
(5)
where I is the current intensity, A V is the potential difference applied, 1 is the length of the domain, and S is the area of each surface where a first-type (Dirichlet) boundary condition is applied. The contribution to the bulk conductivity of the conductance along the surfaces is computed by using the same equation 5 where I is now the total electrical current flowing in or out of the triangulated surface network at the locations of the Dirichlet boundary conditions. 3. S I M U L A T I O N RESULTS" P O R O U S M E D I U M A N D P R I M A R Y D R A I N A G E
As a test material, we selected the UK Permo-Triassic Sandstones (Sherwood S.), which constitute the second largest aquifer system in the UK. This type of porous medium is representative of a fine to medium, poorly consolidated sandstone, with a small fraction of clay and Fe oxide material coating the silicate grains. Consequently, the Sherwood Sandstone is expected to exhibit a non-null solid phase/surface conductivity, suitable for analysis with the model described in Section 2. Samples of the Sherwood Sandstone from South Yorkshire [23] were analyzed for porosity and grain-size distribution. We assumed a log-normal distribution for the grain size distribution and, thus, simulated a gravitationally stable sphere packing that follows the logarithmic mean, standard deviation, and porosity. Simulated packing properties are listed in Table 1. The primary drainage curve was experimentally measured on three remolded Sherwood sandstone samples with a tempe pressure cell apparatus. The non-wetting phase was air
239 Table 1 Properties of the simulated sphere-packing. Par amet er Value U nit grain radius Rg 0.154 + 0.051 mm domain length L 7.974 mm porosity 9 0.39 number of spheres 14501
and the wetting phase was distilled water. The surface tension between air and water, "7 has been assumed equal to 72 x 10 -a N/m. The porous medium described in Table 1 was discretized with increasing resolution of 1283, 2563 and 3843 voxels to simulate the primary drainage curve. Figure 1 compares experimental results with simulated curves at three discretization levels. Results are affected by resolution only at very low water saturation, while the overall shape of the drainage curve is not substantially affected by the discretization level. Analogous discretization effects on simulated p ~ - s ~ curves were observed and discussed in the literature [20]. Our simulations reproduce correctly the features of the experimental data. Particularly, we obtain an excellent agreement for the flat portion of the pC - s ~ curve and for simulated and experimental air entry pressure value, which are related to the uniformity of the medium and to the mean pore throat diameter, respectively. This implies that the simulated packing we used is well representative of the real pore structure.
9
o O <) 9 9 V__
~-250m-: 3:~ / ] E
~aooJ: ~....
9
Simulated data, 1283 voxels Simulated data, 2563 voxels Simulated data, 3843 voxels Experimental data, 1st set Experimental data, 2 nd set
Experi___ment_aldata__,3:dset_j__
lS~ a,IO%
~
.,,
/
~+'~' ~ .o,~+o~,.,~,"
%
o'.a
o'.4
o'.+
Water saturation
o'.a
-'r
Figure 1. Discretization effects on primary drainage curve and comparison with three sets of experimental data.
240 4. S I M U L A T I O N
RESULTS" BULK ELECTRICAL
CONDUCTIVITY
We computed the bulk electrical conductivity of the porous medium for each equilibrium air-water distribution modeled with the morphological simulator. The three-dimensional simulator discussed in Section 2 was set to run on an IBM SP RS/6000 Power4. 4.1. C a l i b r a t i o n of m a t r i x e l e c t r i c a l c o n d u c t i v i t y Two approaches are possible for the description of conductance related to the clay/oxide fraction of the porous medium" 1. surface conductivity is represented by an equivalent conductivity of the whole solid matrix volume, a~; 2. surface conductivity is represented as a proper interfacial property along the solid/water interfaces, cq. We investigated both approaches. In either cases, the unknown parameter needs to be determined by matching the experimental value of saturated bulk conductivity of the sandstone [23], equal to a b - 1.56 x 10 .2 S/re. 4.1.1. solid e q u i v a l e n t c o n d u c t i v i t y a p p r o a c h This calibration was performed on a completely saturated porous medium, discretized with 3843 voxels. We assume aw - 6.8 x 10 -2 S//m, as from experimental measurements by [23] for natural groundwater at the site. Figure 2a shows the results. The resulting bulk solid conductivity is a~ - 1.43 x 10 -3 S/re. As expected, this is a nonzero value, even though significantly lower than water conductivity. In this approach, ~ represents the additional surface conductivity due to mobile ions in the electrical double layer of the clay.
0.019
0.016-
0.0158 ~
0.018
,,,,~' ,,
0.0156 . . . . . . . . . . . . . . . . . . . . . . . .
,,,
0.017
,,'0'"
,
,.-- 0.0154 E
'~'~.016 _
,,,"
,,"'"
"-".,,0.0152 "5 0.015
,,'"
0.015 ~ 0.0148-
0.014
,,,
,,"
,,~'
,,,'" O"
0.0146-
a)
o.o~3o
0'.5
1
,I', 1.5 ~solid (S/m)
2.5
3
x 0 b)
5
6
7
o,s
x 04 '
Figure 2. Calibration of solid phase conductivity for a fully saturated medium with the solid equivalent conductivity approach (a) and the interface conductivity approach (b).
241 Different models have been proposed to take into account the electrical double layer conduction. The effective bulk conductivity of a saturated medium can be defined [10] as
1
~v - T ( B Q v + ~ )
(6)
where Q~ is the volume concentration of clay exchange cations (meq/mL) and B the equivalent conductance of the clay counterions (cm2S/meq), which is related to their mobility. Some works [12] define B Q ~ - (Tclay and compute it assuming conduction pathways through the solid phase. According to this model, our cr~ - cr~lay. The cr~ value we compute is lower than those presented in both papers, indicating a really low surface conductance. For simulations of electrical charge flow under unsaturated conditions, we assumed cr~lay = BQ~ ~ [10]. This states that when water saturation is less than unit, the exchange cations associated with the clay become more concentrated in the remaining pore water.
4.1.2. interface conductivity approach We made also a first attempt at calibrating the surface conductivity according to the second approach, on a small porous medium derived as a 15a voxels subdomain of the digital porous medium discretized with 128 voxels per side. Figure 3 shows the triangular mesh obtained with the marching cubes code that models the interface between solid and aqueous phase.
Figure 3. Triangular mesh used to compute the surface conductivity.
242
Calibration results are shown in Figure 2b. We calculated ai - 9.0 x 10 . 4 S. This value is of the same order of magnitude of the one calculated for a volume equivalent conductivity of the solid phase. Due to the preliminary nature of the results a more specific comparison between the two values may not be meaningful. Particularly investigation of domain size effects on simulated ai deserves further studies. Thus, at this stage, comparison against experimental data has been performed using approach 1. 4.2. F i c t i t i o u s e l e c t r i c a l c o n d u c t i v i t y of g a s e o u s p h a s e An effective value for electrical conductivity of the gaseous phase (air) has to be determined for numerical reasons, since a value equal to zero would possibly result in a singular matrix of the finite difference system of equations. We adopted a~ - 6.8 • 10 -2 S / m and a~ - 1.43 • 10 -3 S/m and calibrated the value of gaseous phase electrical conductivity on a porous medium discretized with 3843 voxels with s~ - 0.13 (a high air content). The gas conductivity, physically zero, must be chosen as large as possible in order not to generate a matrix too ill-conditioned, but must be small enough so that any smaller gas conductivity would effectively yield the same medium bulk conductivity (asymptotic behavior) (see Figure 4). The asymptote is reached for O"a nearly equal to 1.0 • 10 -7 S/m, typical of an insulator. Note that the computational time increases with decreasing gaseous phase conductivity~ as the system matrix condition number increases, and forces to the iterative solver to converge more slowly.
x 10 -4 ,
i
9.54 9.52
....-
9.5
.-
9.48 E (/)
9.46
,,,,,,,'~
9.44 ,3
9.42 9.4 9.38 9.36 10 -7
lO-(s/m )
10 -5
(~air
Figure 4. Calibration of gaseous phase conductivity.
4.3. C o m p a r i s o n w i t h e x p e r i m e n t a l d a t a The results of electrical field simulations on the 3843 voxels porous medium are shown in Figure 5 as a function of water saturation~ together with experimental measurement
243 results [23]. Comparison between simulated and experimental bulk conductivities show an excellent agreement at high and intermediate s ~ values.
450 400
I
- O " Simulated data 9 Experimental data
350 300 E
~250 -5 .o (:2.
200 150 100 1
0.4 0.6 Water saturation
i
0.8
1
Figure 5. Simulated bulk conductivities and experimental data from [23] as a function of water saturation.
5. C O N C L U S I O N S
A pore-scale method based on digital representation of a porous medium was developed and applied to predict the bulk electrical conductivity of unsaturated porous media. The model was used to simulate the DC resistivity of semi-consolidated Triassic sandstone from the UK (Sherwood sandstone). Comparison with experimental data show a very good agreement with model results. The results presented in this paper show the accuracy and reliability of the proposed method for pore-scale modeling and computation of bulk conductivities of a porous medium. This method can be a powerful tool in the investigation of constitutive electrical law from a pore scale point of view. Further work will focus on testing the surface conductivity model and extend the model to account for induced polarization phenomena. REFERENCES
1. O. Mazac, W. E. Kelly, I. Landa Surface geoelectrics for groundwater pollution and protection studies Journal of Hydrology, vol.93, no.3-4, pages.277-294, 1987. 2. R.K. Frohlich and C. D. Parke. The electrical resistivity of the vadose zone - field survey. Ground Water, 27(4):524-530, 1989. 3. W. D. Daily, A. L. Ramirez, D. J. LaBrecque, and J. Nitao. Electrical resistivity tomography of vadose water movement. Water Resour. Res., 28:1429-1442, 1992.
244 4.
R. J. Kalinski and W. E. Kelly. Estimating water content of soils from electrical resistivity. Geotechnical testing journal, pages 323-329, 1993. 5. A. Ramirez, W. Daily, A. Binley, D. LaBrecque, and D. Roelant. Detection of leaks in underground storage tanks using electrical resistence methods. J. Env. and Eng. Geophysics, 1(3)'189-203, 1996. 6. A. Binley, G. Cassiani, R. Middleton, and P. Winship. Vadose zone model parameterisation using cross-borehole radar and resistivity imaging. J. Hydrol, 267(3-4)" 147-159, 2002. 7. D. Michot, Y. Benderitter, A. Dorigny, B. Nicoullaud, D. King, and A. Tabbagh. Spatial and temporal monitoring of soil water content with an irrigated corn crop cover using surface electrical resistivity tomography. Water Resources Research, 39
(5), 2003. 8. 9. 10. 11.
12. 13.
14. 15.
16.
17. 18. 19. 20. 21. 22.
G. E. Archie. The electrical resistivity log as an aid in determining some reservoir characteristics. Trans. Am. Inst. Min. Metall. Pet. Eng., 146:54-62, 1942. A.E. Bussian. Electrical conductance in a porous medium. Geophysics, 48:1258-1268, 1983. M. H. Waxman and L. J. M. Smits. Electrical conductivities in oil-bearing shaly sands. Soc. Pet. Eng. J., 243:107-122, 1968. C. Clavier, G. Coates, and J. Dumanoir. The theoretical and experimental bases fro the dual-water model for the interpretation of shaly sands. Soc. Pet. Eng. J., 24(2)" 153-169, 1984. P. W. J. Glover, M. J. Hole, and J. Pous. A modified archie's law for two conducting phases. Earth and planetary sciences letters, 180:369-383, 2000. A. Revil, L. M. Cathles III, and S. Losh. Electrical conductivity in shaly sands with geophysical applications. Journal of geophysical research, 103(B10)'23,925-23,936, 1998. A. Yang, C. T. Miller, and L. D. Turcoliver. Simulation of correlated and uncorrelated packing of random size spheres. Physical Review E, 53(2)'1516-1524, 1996. A. L. Endres and R. Knight. The effects of pore-scale fluid distribution on the physical properties of partially saturated tight sandstones. J. Appl. Phys., 69(2)'1091-1098, 1991. R. J. Suman and R. J. Knight. Effects of pore structure and wettability on electrical resistivity of partially saturated rocks - a network study. Geophysics, 62(4)'1151-1162, 1997. D. Zhou, S. Arbabi, and E. H. Stenby. A percolation study of wettability effect on the electrical properties of reservoir rocks. Transport in porous media, 29:85-98, 1997. C. D. Tsakiroglou and M. Fleury. Pore network analysis of resistivity index for waterwet porous media. Transport in porous media, 35:89-128, 1999. H. N. Man and X. D. Jing. Pore network modelling of electrical resistivity and capillary pressure characteristics. Transport in Porous Media, 41:263-286, 2000. M. Hilpert and C. T. Miller. Pore-morphology-based simulation of drainage in totally wetting porous media. Advances in Water Resources, 24:243-255, 2001. Lorensen, W. E. and Cline, H. E. Marching cubes: A high resolution 3D surface construction algorithm. Computer Graphics, 21:163-169, 1987. E. Dalla, M. Hilpert, and C. T. Miller. Computation of the interracial area for multi-
245 fluid porous medium systems. Journal of Contaminant Hydrology, 56(1-2)'25-48, 2002. 23. A. Binley, P. Winship, R. Middleton, M. Pokar, and J. West. Observation of seasonal dynamics in the vadose zone using borehole radar and resistivity. In Proceedings
of the symposium on the application of geophysics to environmental and engineering problems, SAGEEP, march 4-7 2001.