Soil is a complex, potentially heterogeneous, lossy, and dispersive medium. Modeling the propagation and scattering of electromagnetic (EM) waves in soil is, hence, more challenging than in air or in other less complex media. This chapter will explain fundamentals of the numerical modeling of EM wave propagation and scattering in soil through solving Maxwell’s equations using a finite difference time domain (FDTD) method. The chapter will explain how: (i) the lossy and dispersive soil medium (in both dry and fully water-saturated conditions), (ii) a fourth phase (anomaly), (iii) two different types of transmitting antennae (a monopole and a dipole), and (iv) required absorbing boundary conditions can numerically be modeled. This is described through two examples that simulate the detection of DNAPL (dense nonaqueous-phase liquid) contamination in soil using Cross-well radar (CWR). CWR —otherwise known as cross-borehole GPR (ground penetrating radar)—modality was selected to eliminate the need for simulation of the roughness of the soil-air interface. The two examples demonstrate the scattering effect of a dielectric anomaly (representing a DNAPL pool) on the EM wave propagation through soil.
The objective behind selecting these two examples is twofold: (i) explanation of the details and challenges of numerical modeling of EM wave propagation and scattering through soil for an actual problem (in this case, DNAPL detection), and (ii) demonstration of the feasibility of using EM waves for this actual detection problem.
In addition, the results of an experimental simulation of Example 1 (the case with the monopole antenna) will be analyzed and discussed. The results of the corresponding numerically simulated example will then be compared and validated against the above-mentioned experimental results. A conclusion section will close the chapter.
Before explaining the numerical modeling and its challenges, some background about DNAPLs and their detection technologies, including the CWR method, is explained in the following section.
Cross-well radar (CWR), otherwise known as cross-borehole ground-penetrating radar (cross-borehole GPR) is a minimally invasive method that uses high frequency electromagnetic (EM) waves transmitted and received by antennae in the subsurface to image objects of contrasting dielectric properties. In order to assess the feasibility of using CWR for detection of dense nonaqueous-phase liquid (DNAPL) pools in deeper layers, this technique should be numerically simulated and evaluated. CWR is more appropriate for deep investigations. For near-surface sensing such as bridge-deck health monitoring, methods such as GPR are more appropriate (Belli et al., 2009 and 2009a). GPR studies require addressing the scattering due to the rough soil-air interface in the forward model as well as through the inversion process (Firoozabadi et al., 2007).This chapter describes a numerical modeling approach to Maxwell's equations using a finite difference time domain (FDTD) solution with both monopole and dipole antennae to simulate the scattering and propagation of EM waves in DNAPL-contaminated media. An FDTD code, originally developed for detection of mines using two-dimensional (2D) surface-reflection ground-penetrating radar (GPR) through non-dispersive media, was revised and upgraded for three-dimensional (3D) cross-borehole wave propagation in heterogeneous soils. The three dimensional FDTD code was enhanced to accommodate dispersive media and was used to model EM wave propagation in heterogeneous soils. This chapter describes the effect of the radiation patterns of two different antennae on propagation and scattering of EM waves through the soil subsurface and its potential for detection of DNAPL pools. In order to evaluate the feasibility of using the CWR method to detect DNAPL pools, illustrative examples with and without the presence of the DNAPL pool were analyzed. The results show considerable diagnostic potential to detect contaminated zones with DNAPLs using EM waves through CWR.
DNAPLs are separate-phase hydrocarbon liquids denser than water, such as chlorinated solvents (tetrachloroethene (PCE) and trichloroethylene (TCE)), wood preservatives, coal tar wastes, and pesticides. DNAPLs may not usually be found as a free phase in soil cores or accumulated in monitoring wells. Based on this lack of observable evidence of DNAPL pools, investigators may conclude that no DNAPL is present, when it may be present in substantial quantities at residual saturation as large as 70% to 80% of the total porosity (ITRC, 2000).
There are varieties of invasive techniques to detect DNAPLs such as direct push probe techniques (e.g., direct soil sampling or indirect sampling such as negative ion sensor) and use of in situ tracers (e.g., PITT, or partitioning interwell tracer test), excavating test pits, and groundwater profiling. Most invasive techniques just provide point-sources of information and may help the DNAPL pool to spread through the substantial number of required boreholes and excavated trenches, while noninvasive geophysical techniques noninvasively or minimally invasively detect DNAPLs. They can use different types of waves (e.g., electromagnetic, acoustic, etc.) to indirectly reconstruct images to characterize or detect anomalies as DNAPLs within heterogeneous media. These methods are minimally invasive techniques that discriminate the contrast between local physical properties of the background and target to produce images of the subsurface.
Cross-well P-wave transmission at 90 kHz was used in a sand pack before and after introducing NAPLs by Geller et al. (2000). The results indicated that small NAPL saturations may be more easily detected with amplitude than with travel time data, but the relationships between the amplitude changes and NAPL saturation may be more complex than those for velocity.
Cross-hole complex resistivity was also applied to a contaminated vadose zone by Grimm and Olhoeft (2004) to predict the general distribution of DNAPLs at parts per thousand concentrations, specifically widespread near-surface contamination and in the vadose zone immediately underneath the source.
Complex resistivity (CR) is a technique that so far has shown a potential in detecting DNAPL pools (Blackhawk Geoservices Inc., 2008).
Smith-Rose (1993& 1935) is among the earliest scientists who studied the dielectric permittivity of different soils. Some physical properties of a typical DNAPL (i.e., PCE) are compared with those of different soils, water, and air in Table 1.
As seen, dielectric permittivity and effective (= DC + alternating field) electrical conductivity of PCE are lower than those of water and relatively similar to air or some dry soils. Thus, a PCE-contaminated region of the saturated zone should provide a contrast with a similar magnitude of the unsaturated-saturated zone interface (water table). The values of dielectric permittivity for water and soil are dispersive (depend on frequency), but most DNAPLs are nonpolar molecules with minimal frequency dependence and therefore with almost no dispersive characteristics. In frequencies greater than 10 MHz and less than 1 GHz, the dielectric permittivity is controlled by the polarization of individual water molecules and is therefore dependent on moisture content (Binley et al., 2001). Dielectric permittivity in the lower MHz range of frequency (< 10 MHz) depends on particle shape and mineralogy, electrolyte type and concentration, particle orientation, and soil electrolyte interaction (Rinaldi & Francisca, 1999). In the upper MHz range of frequency (10 to 1000 MHz), the real part of dielectric permittivity is affected by the polar contribution of bound and free water molecules (Hoekstra & Doyle, 1971; Hoekstra & Delaney, 1974; Selig & Mansukhani, 1975; Dobson et al., 1985; Hallikainen et al., 1985; and Arulanandan, 1964). Sachs and Spiegler (1964) presented a model for dielectric permittivity of soil mixtures at different frequencies. This model is based on an equivalent circuit for conductive particles in electrolytes and fitted in the dielectric plane by adjustment of three parameters. Arulanandan and Smith (1973) explained some aspects of dielectric dispersion based on the Sachs and Spiegler (1964) model in a frequency range from 106 to 108 Hz. Another model that explains the dispersion of the soil relative dielectric permittivity, otherwise known as dielectric constant, was presented by Thevanayagham (1995), which considers the effect of particle orientation relative to the electric field. Weedon and Rappaport (1997) and Rappaport et al. (1999) presented a model to predict the dispersive nature of the dielectric permittivity and electrical conductivity of soils, using a rational Z-transform approximation function of the conductivity. This technique uses the values of dielectric permittivity and electrical conductivity of the bulk mixture, and does not deal with the components and their volumetric content. It is basically a technique to convert the bulk values for the matrix to a form applicable to the FDTD method to take the dispersive nature of the resultant mixture into account.
|Material||Density ρ (g/cm3)||Real Part of Dielectric Constant* (εr)||**** Effective Electrical Conductivity σe (1/Ω.m)|
|f = 100 MHz||f = 1.5 GHz||f = 100 MHz||f = 1.5 GHz|
|PCE||1.62||2.28||2.28||2.63 × 10-6||1.52 × 10-4|
|Water**||1.00||78||77.14||2.2 × 10-3||0.509|
|Dry Sand||1.30||2.55||2.55||1.8 × 10-4||2.1 × 10-3|
Dielectric constant is also known as relative dielectric permittivity** Pure water *** 20% moist Bosnian clay loam**** Effective electrical conductivity (σe) = Static or DC conductivity (σs) + Alternating field conductivity (σa)
The loss (attenuation) in soil is a function of a variety of factors such as soil type and mineralogy, moisture content, electrical conductivity, and frequency. Higher frequency EM waves attenuate faster. In other words, soil acts as a low-pass filter. However, higher frequency EM waves provide higher image resolution. There is a tradeoff between the penetration depth and image resolution. Therefore, feasibility of the CWR to detect an average size DNAPL pool using proper frequency range is the first step to detection, which is the main goal of this chapter.
Soil, water, and DNAPL mixed at DNAPL saturations of 0% to 50% (i.e., water saturation of 50% to 100% below water table; i.e., no air) have been evaluated by Ajo-Franklin et al. (2004) to have a dielectric constant of 9 (50% water saturation soil with 50% DNAPL saturation) to 24 (fully (100%) water-saturated soil with 0% DNAPL saturation). These values are different from the dielectric constant of pure DNAPL (ε’ ≈ 2.3) and vary based on porosity and degree of DNAPL saturation as the result of mixture.
There are some expected limitations to using CWR to detect DNAPLs in soil and groundwater such as: (a) use of CWR to detect DNAPLS may not perform as well in dry or partially saturated zones as in water saturated zones, due to the weaker contrast in electromagnetic properties of DNAPLs and dry soils compared to the one between DNAPLs and water saturated soils, (b) the concentration (or contaminant saturation) must be fairly high, which makes the detection of DNAPLs using CWR often limited to identifying DNAPL pool sources (high DNAPL saturation) and not plumes, and (c) the spacing between the wells strongly influences the effectiveness of CWR (as the separation between transmitting and receiving antennae increases, the radar wave amplitude attenuates, which creates greater difficulty in distinguishing the wave from background noises) (ITRC, 2000).
Therefore, in this chapter a simple case of DNAPL dominantly replacing water in the pores (DNAPL saturation of 80% to 100%) is simulated numerically. Based on Bruggeman-Hanai-Sen (BHS) model (Sen et al., 1981), due to low dielectric constant of soil grains (≈ 5, (Ajo-Franklin et al., 2004)), the bulk dielectric constant of the soil mixture within the DNAPL pool (close to 100% DNAPL saturation) values between 2.3 and 5. Hence, for comparison, acrylic was selected as the dielectric object of ’ ≈ 2.6 (Weast, 1974) since it has been used in another experimental work by the authors (Farid et al., 2006). This strong contrast in dielectric constant between the dielectric constant of the acrylic (≈ 2.6) and water-saturated soils (≈ 15 - 25) makes detecting the scatterer in the saturated soil using radar-based geophysical methods feasible.
The U.S. Geological Survey used 500-MHz surface GPR among other methods to monitor the location and migration of the subsequent plume (Sneddon et al., 2000).
Bradford and Wu (2007) tested 3D multi-fold GPR on a small scale controlled DNAPL release to detect contaminated soil with demonstrated success.
While GPR is the least invasive radar-based method, it is not a practical one for deep investigations. Cross-Well Radar (CWR), otherwise known as cross-borehole GPR, may be a more effective method for deep contaminant detection. CWR is a method that uses EM wave transmission measurement across borehole antennae as opposed to reflection measurements used in surface-reflection imaging methods, such as GPR. CWR overcomes some of depth limitations of surface-reflection imaging GPR. CWR uses antennae that are lowered into sampling wells with cables. Radar waves are emitted from a transmitting antenna in one well and propagate through the ground to a receiving antenna in a second well. The subsurface geology and pore fluids scatter, refract, or reflect the waves. DNAPL pools also provide dielectric permittivity and electrical conductivity contrasts within relative moist soils, which in turn, cause refractions or reflections on EM waves. These scatterings and reflections make DNAPLs amenable to detection by radar. Other than the magnitude of the scattering field, phase, and hence, travel time can also be used for other techniques such as travel time tomography to reconstruct images of the background medium and anomalies within the medium.
Success of radar-based methods as well as any other detection technique depends on the strength of the contrast in properties of the background and scatterer (e.g., DNAPLs, moisture). Dielectric permittivity and electrical conductivity are the main parameters controlling such differences for radar-based techniques.
CWR relies on a one-way travel time whereas the surface-reflection imaging GPR relies on the 2-way travel time. Therefore, the effective antenna separation of CWR can be twice as far as the penetration depth of surface-reflection GPR, theoretically. GPR antennae are designed to be relatively impedance-matched to the ground, and the reflection coefficient at the air-soil interface is close to unity leading to the potential for strong internal reflection. However, there is still the disadvantage of reflections at the air-soil interface, which reduces the amount of received waves reflected by the underground objects. The performance of detection using CWR is highly dependent on seasonal conditions, in particular soil moisture (Daniels et al., 1992). Soil mineral composition and physicochemical properties also influence radar signals. For example, clays tend to attenuate these signals and therefore may limit the skin depth and effectiveness of the method (Anderson & Peltola, 1996). Considering these limiting factors, feasibility of the technique for DNAPL detection should be evaluated by simulating implementation of ultra-wideband CWR. Further understanding of the behavior of 3D EM waves in soil is necessary for the future implementation and feasibility evaluation of CWR for detecting DNAPLs. To achieve this goal, the detection technique using monopole and dipole antennae is modeled numerically in the time domain using an FDTD technique. A critical issue addressed in this research, by direct simulation of the antennae, is the wave interaction at the antenna/soil interface. The technique models the transmitting antenna but measures the electrical field at all grid points on any desired cross-sectional or depth slice (instead of modeling receiving antennae). This is equivalent to having hundreds of receiving antennae in soil, which is not practical. In practical techniques such as cross-well tomography, few antennae are installed and used alternatively as transmitters or receivers to collect data in a multiple-depth, multiple-location manner, and the outcome is used for inversion and image processing techniques (the authors are working on different aspects of cross-well tomography and results will be published in the future). Obviously, the simulation technique presented in this chapter is used only for preliminary evaluation of the feasibility of detecting DNAPL pools using CWR. The outcome can also be used as a forward modeling for future inversion techniques.
Besides, another goal of this chapter is to evaluate the feasibility of DNAPL detection and to investigate if there are points with strong signature by the object on the scattered field, not to discuss inverse modeling. Therefore, receiving antennae are not modeled. The results model no noise. Nevertheless, comparison between the values of the scattered field (by the target) with those of the incident field will provide helpful information that can lead to the critical noise/signal ratio, even without applying any noise. This is discussed in more detail in the following sections. Experimental works are being conducted to evaluate typical levels of noise, compare them to the levels derived from this technique, then apply the noise to the simulation and reevaluate the feasibility (Farid et al., 2006).
3. FTDT modeling of maxwell’s equations
Finite Difference Time Domain (FDTD) is a relatively powerful and very popular method because of its simplicity and despite its instability problems. Therefore, a 3D FDTD code was developed and used to simulate the antenna performance, coupling, and interactions in soils. The FDTD technique was originally introduced by Yee (1966) and is based on time and spatial discretization of Maxwell’s equations to obtain solutions for the EM field in the time domain. The technique is numerically implemented by continuously sampling the electromagnetic field over the wave propagation medium discretized into a grid. The differential form of Maxwell’s equations in the time domain is as follows (Grant & Philips, 1990; Sheriff, 1989; and Balanis, 1989).
along with the constitutive relations,
where E is the electric field density, H is the magnetic field intensity, D is the electric displacement or electric flux density, B is the magnetic flux density, is the impressed (source) electric current density, is the impressed (source) magnetic conductive current density, e is the electrical conductivity, m is the magnetic resistivity, ε is the dielectric permittivity, is the magnetic permeability, ρe is the electric charge density, and ρm is the (monopole) magnetic charge density. All of the field parameters, E, H, B, D, , and are assumed to be functions of position and time t, while material parameters , , and are functions of position and may be dispersive functions of frequency (and thus, through the Fourier Transform, also of time). Note that the usual multiplicative constitutive relations are written as convolutions for dispersive soil media. Equations (1) through (4) can be also written in integral forms. These four partial differential equations are discretized in a 3D rectangular grid, using a central-difference approximation of two consecutive values for the field components in both space and time. Due to the nature of electric and magnetic fields (e.g., coil loops), these two cannot be located at the same spatial location. In the discretization process, electric and magnetic field components are assigned to the edges of complementary interlocking cubical meshes. Thus, half-space indices are introduced for different field components. Besides, half time steps are introduced to perform the finite difference computation of E based on H and vice versa. Spatial and time steps are respectively represented by the lower indices (i, j, k) and the upper index (n). Some details of the simulation process are explained in the following. Fig. 1 helps to better understand this staggered time and space grid. As seen in the superscripts in the figure, there is ½ time step difference between E and H. Due to the central difference approximation technique in time, H is present at both the (n + ½)th and (n - ½)th time steps in the equations along with E only at the nth time step.
3.1. Advantages and limitations
The most advantageous characteristic of the FDTD method is its simplicity. Solving Maxwell’s equations using FDTD is a simple iterative procedure. However, there are several limitations to the numerical implementation of the FDTD technique through the above-mentioned difference equations. These limitations have been and are subject of research. Some of these challenging limitations impose restrictions on the grid size and time step increments that affect stability and accuracy of the method (Kunz & Luebbers, 1993).
One fundamental restriction requires the longest side of the grid cell to be much shorter than the shortest wavelength of the wave within the cell. A very common restriction assumed in practice is / 10, where is the shortest appreciable wavelength in the excitation signal. This constraint is imposed by both sampling limitation and grid dispersion errors (Kosmas, 2006).
The second restriction develops from the scale and geometry of the problem. The method uses a uniform grid to model small antennae along with a large soil medium with or without contrasting objects. Accordingly, the geometry imposes challenging limitations, especially in computation cost. For example, in the case of a DNAPL contaminated site, dimensions of the background in horizontal directions may be of the order of tens of meters by tens of meters in the XY-plane, while the antenna thickness may be of the order of millimeters to centimeters. This requires solving a large uniform grid with a small grid-cell size in the X and Y directions. At the same time, a small vertical (Z-direction) cell size is required to account for thin DNAPL pools. One way to solve this problem is by using a non-uniform grid, which adds more difficulties to the task of satisfying stability conditions (e.g., Courant condition).
The third limitation is the time-step restriction required to satisfy the Courant condition. The Courant condition restricts time increments to a range in which waves do not travel too far in each time increment, or:
where N is the grid dimension, xm is the grid cell size in the mth direction, and Vp,max is the maximum wave phase velocity within the model. In a 3D case, Equation (9) will be written in the following form:
In the case of N-dimensional isotropic cells, Equation (9) can be simplified to (e.g., in a three dimensional case, N = 3, and x = y = z = ), where the left hand side is called the Courant number. Different cases have demonstrated that using smaller values of t does not necessarily improve the results. However, smaller values for the Courant number may sometimes yield satisfactory results (Kosmas, 2006).
4. Modeling procedure
4.1. Modeling dispersive nature of soil
Soil is a complex medium to model because of its heterogeneous, lossy, and dispersive nature. Accurate computation of the behavior of soil over a wideband of frequency requires either several individual frequency domain calculations or a robust deconvolution of E(t) from D(t) in the time domain (Rappaport et al., 1999). The convolutional approach to modeling dispersion in soil approximates the frequency domain dispersive complex dielectric constant with rational functions of jω and multiplies the constitutive relations by the denominator. Then, the results are inverse-Fourier-transformed. This model is called the Debye or Lorentz model (Kashiwa & Fukai, 1990; and Gandhi, 1993). Weedon and Rappaport (1997) and Rappaport and Winton (1997) simplified the problem by modeling the conductivity as a simple rational function of Z-transform. Modeling conductivity in terms of the Z-transform variable (Z-1 readily transforms to “time delays”) simplifies the process of conversion of the generalized dispersive Ohm’s Law (J(Z) = σ(Z) E(Z)) to the time domain. To keep the simplicity of the FDTD method, Z-transform function of conductivity in the time domain is given as (Weedon & Rappaport, 1997; and Rappaport et al., 1999):
Using this function, both real and imaginary components of conductivity depend on the sampling interval Δt and the coefficients of the rational Z-transform function (a1, b0, b1, and b2). The real part of conductivity corresponds to measured conductivity, and the imaginary part corresponds to real dielectric permittivity. Ampere’s law can be written in the Z-domain as follows.
4.2. Absorbing boundary conditions
As in simulation of any diffusion application (e.g., heat transfer or water flow through porous media), wave propagation through infinite media is practically modeled in a finite grid. To model the infinity of the flow or wave propagation through the required infinite media by the finite number of grid cells, appropriate boundary conditions are required. The lattice termination absorbing boundary condition used in the FDTD code was based on the second-order one-way wave propagation equation of Mur (1981) and adjusted for lossy soils by Talbot and Rappaport (2000) using the rational function approximation of Equation (11).
4.3. Antenna and excitation
There are some factors that categorize the excitation, such as the source geometry, excitation signal type, and hard or soft source of propagation.
Hard and Soft Sources: The code uses soft sources for all different antenna types. A hard source specifies the total field at the excitation point. These types of sources are usually avoided as they cause undesirable reflections at the source point. Alternatively, a soft source specifies the additional field supplied above the existing background field at the source point. The soft source can be specified at points on a radiating aperture, while hard sources are used for current sources flowing along metal structures. In 2D cases, soft sources propagate well, but extensive testing of 3D FDTD cases has shown difficulties involved in propagation from a soft source unless the physical body of the antennae is modeled. This paper uses soft sources for both monopole and dipole antenna cases of the following sections.
Excitation Signal: A cosine-modulated Gaussian signal (Equation (13)) was used to excite both antennae because of its simplicity and frequency content
where E0 is the amplitude of the pulse, t is any arbitrary time instant, and t0 is the time instant corresponding to the peak of the Gaussian. Ideally, a Gaussian pulse has infinite duration. However, in numerical implementation it ought to be truncated. Thus, W is chosen as the duration of the pulse. Parameters t0 and W can be chosen, but not absolutely arbitrarily. In this chapter, W was chosen based on the bandwidth used in a scaled experiment by the authors (Farid et al., 2006). Inappropriate choices of these parameters may introduce instability and noise into the FDTD computation. Starting and truncation of the Gaussian signal should be sufficiently smooth (i.e., the field value at the truncation time should be sufficiently small) to alleviate introduction of frequencies with spectrum levels above the one that a single or double precision calculation can tolerate. Large spectrum energy of the incident field may cause instability and noise, which in turn, results in computation corruption. This is true for any type of pulse. For the Gaussian type pulse, experience has shown that a choice of t0 ≥ 4W ensures smooth truncation, and in turn, good performance (Kosmas, 2006).
4.4. Geometry of the source
To model the DNAPL detection and evaluate its feasibility, two different antenna types were modeled, a monopole and a dipole. Some fundamental parameters and facts about these simulations are explained in the following.
The monopole antenna is simply a coaxial cable, schematically shown in Fig. 2(a) with its shield removed at a specific length from the tip (referred to as extended dielectric length). In a simple monopole antenna, the most efficient radiation occurs when the exposed dielectric and center conductor length is ¼ of the excitation wavelength. For the range of frequency of interest of this and other ongoing and further research (400-2200 MHz), and based on the dielectric permittivity of the dielectric part of the antennae (≈ 2.1), this length should be between 2.35 cm and 12.94 cm. Different types of antennae have different radiation patterns. The radiation pattern is a visual way of representing how an antenna distributes energy into the surrounding space. Different methods (parameters) can be used to demonstrate the wave propagation through space. One popular method is to use power as the key parameter. This method visualizes how and in which direction power is concentrated by the antenna to propagate away. This is possible by drawing contours of equal power. Other parameters, such as field magnitude or phase, can be used for visualization purposes. Contours of these parameters are drawn to visualize radiation patterns. A monopole antenna radiates waves in all different directions.
A dipole antenna consists of two straight metallic or dielectric parts connected at the center to a feed line. A typical dipole antenna is shown in Fig. 2(b). This antenna type constitutes the main RF radiating and receiving element in various sophisticated types. Dipoles are inherently balanced antennae due to bilateral symmetry. The best dipole length is usually ½ wavelengths.
5. Evaluating feasibility of DNAPL detection
As mentioned before, two different antennae are simulated in the code to model the problem and evaluate the feasibility of the use of the CWR method to detect dielectric objects in soil. The first case models a small-scale modeled monopole antenna as a soft source of the wave within a fully water-saturated sandy soil medium (degree of water-saturation = Sr = 1) and the gravimetric moisture content = w = 17%). The second case is the simulation of a dipole antenna in a larger size medium of the same fully saturated sandy soil (Sr = 1 and w = 17%). This value of moisture content is selected to model a soil medium similar to the existing fully saturated sandy soil in a pilot-scale facility constructed at Northeastern University (referred to as SoilBED) for subsurface sensing and imaging experimentation. This 17% moisture content is the reference moisture content for much further experimental research by the authors (for more information refer to Farid et al., 2006) in the SoilBED facility. Hydrological modeling of DNAPLs and the dielectric properties converted through a generic petrophysical relationship in the simulation is desired but does not fit in the extent of this chapter. Details of the two above-mentioned simulated cases follow.
5.1. Monopole antenna case (A)
The pilot-scale simulation of the SoilBED facility (Farid et al., 2006) is explained in this section. In the first case, a 5 mm-thick monopole antenna is modeled within a fully saturated sandy soil background. The size of the medium under study was selected to satisfy limitations of the FDTD code as well as the experiment. Table 2(a) summarizes details about the geometry and grid size of the soil medium. The simulation is driven by a cosine modulated Gaussian time pulse at a reasonably high frequency (1.5 GHz). To accommodate the simulation of the dispersive soil and stability of the FDTD code at this frequency, the time increment Δt = 2 psec was used. The dispersive properties of the soil for this choice of center frequency and time-step are modeled with the Z-transform function coefficient set, a1, b0, b1, and b2, given in Table 2(b).
To relate the results to the field site, the model can be scaled up in size while scaling down the frequency. To evaluate the feasibility of the DNAPL detection method using monopole antennae, wave propagation through the background soil and scattered EM wave propagation by a DNAPL pool were modeled and analyzed. The geometry details of the monopole transmitting antenna modeled in this case are tabulated in Table 2(c). The drive signal excites the top of the simulated coaxial cable feeding the monopole antenna in a conventional radial field pattern. The electric field components on all grid points of different cross-sectional and depth slices of the medium were computed and then visualized using MATLAB.
First, the background medium was analyzed. Then, a rectangular acrylic plate as a representative of a DNAPL pool was modeled within the soil medium. Fig. 3 schematically shows the simulated geometry (the monopole antenna, the DNAPL pool, and the soil medium). Details of the geometry of the DNAPL pool scatterer are listed in Table 2(d).
|Simulated grid||149 × 149 × 29|
|Grid cell size||0.2 cm × 0.2 cm × 1 cm|
|Entire grid size||29.6 cm × 29.6 cm × 28 cm|
|Soil thickness||21 cm|
|Air thickness||7 cm|
Due to solving the problem at Δt = 2 psec, the FDTD code is very sensitive, and all digits are necessary to satisfy the stability conditions
|Antenna depth||120 mm|
|Perfectly conducting core wire thickness||1 mm|
|Extended dielectric length||20 mm|
|Extended dielectric thickness||3 mm|
|Perfectly conducting outer conductor (shield) thickness||3 mm|
|Gaussian width||0.667 nsec|
|Gaussian peak||5 nsec|
The dielectric constant and effective electrical conductivity of the extended dielectric of the antenna are respectively assumed to be 2.1 and zero (Ω-1).
|DNAPL Pool Geometry||Size|
|Horizontal cross-section||3 cm × 3 cm × 1 cm|
|Clear separation from the antenna||3.8 cm|
|Coordinate of the pool center*||6 cm, 0 cm, -2 cm|
With respect to the center of the grid
To evaluate the wave propagation, the following observation slices were selected. Different components of electric field were computed and visualized on these slices.
A cross-sectional (horizontal: XY-plane) slice, cutting through the antenna and DNAPL pool at the depth of 9 cm. Z and X components of the electric field (Ez and Ex) are shown on this slice in Fig. 4.
Up to this point, only the three vector components of the electric field were visualized. Now, the power is depicted. The intensity of a rapidly varying field is often displayed on a dB scale, enabling the visualization of small amplitude levels. This scale is given by 20 log10|E / Emax|. It is important to note that on the selected depth slice, Ey equals 0, and hence E = Ex+Ez. In addition, since the time domain signals are all purely real, but may have positive or negative values, the dB scale is artificially augmented with positive values to indicate negative field values and better display the oscillating nature of the rapidly decaying wave. The sign of corresponding Ez governs the sign of the dB value. It should be stressed that 0 dB is the maximum field intensity, and positive dB values correspond to weaker signals with the opposite sign.
A depth (vertical: XZ-plane) slice, passing through the antenna and DNAPL pool. This slice (XZ-plane) was chosen because the YZ-plane does not intersect the DNAPL pool. Due to symmetry, Ey is zero on this XZ slice, and hence |E| can be computed by only Ex and Ez. Results are shown in Fig. 5.
This case was initially analyzed without DNAPL contamination (incident field or background) and then with the DNAPL pool (total field). The scattered field by the DNAPL pool target can be computed by subtracting the two previous fields. Three figures are shown for each slice and for each electric field component and for: (i) “incident” (i.e., background, no target), (ii) “total” = background + DNAPL pool target as the scatterer, and (iii) “scattered” (i.e., signature of the target). All results shown in Fig. 4 are captured at t = 3.6 nsec. As seen, incident results of Figs. 4(a) and 4(d) are symmetric, while the total field results shown in Figs. 4(b) and 5(e) are not symmetric. The resulting scattered field information shown in Figs. 4(c) and 4(f) is asymmetric as well.
The incident, total, and scattered (target signature) fields are shown in Fig. 5. The monopole antenna was modeled as a Z-polarized antenna. Therefore, the Z-component of the electric field is the major component, but the scattered field by the DNAPL pool is also readily visible on the X and Y component plots. Since Ez dominates and the scattered field is visible on the Z-component (Fig. 5(c)), the scattered field shown on the dB plot will be clear as well. Further studies (that do not fit in this chapter) show weaker scattered Z-component in dry sandy soils. Different components can be experimentally measured using a receiving antenna with a different polarization (e.g., an X or Y polarized antenna, which is simply a monopole placed horizontally) than the Z-polarized (vertical) transmitting antenna. The scattered field is comparable to the incident field in this case. This potential can also be demonstrated in a different form as shown in Fig. 6.
This figure shows that there is a considerable magnitude and travel time difference between the total and incident fields received at a receiver located right above the DNAPL pool. The strong magnitude difference (more than 100%) and time difference (around 100 psec) between the two signals illustrate the potential of the cross-borehole GPR method to detect DNAPL pools. The early arrival of the total field is caused by the increase in the velocity of EM waves through the DNAPL pool due to its lower dielectric permittivity compared to the saturated soil. The increase in the magnitude of the total field is, on the other hand, caused by lower loss through the DNAPL pool due to its lower electrical conductivity. This illustrates a great potential for DNAPL detection using CWR in saturated soils, if the thickness and size of the pool is a reasonable fraction of the wavelength.
5.2. Dipole antenna case (B)
The above-mentioned small size monopole case can be scaled up to a more realistic size contaminated site. However, scaling up the results may cause some problems that do not allow a simple and direct generalization from small numerical models to real size contaminated sites. For example, in a non-dispersive medium, linear enlargement of the size can be simply interpreted to a linear increase in the wavelength and decrease in the frequency. However, in a dispersive medium, any change in the frequency causes variations in the dielectric properties of the medium. This change in the dielectric constant causes variations in the wave velocity, which in turn adds nonlinearity to the scaling process from the simulated medium up to the real size.
Therefore, to evaluate the scaling issues in a dispersive medium and study the effect of different radiation patterns of different antennae, another case with a more realistic size of soil medium surrounding a dipole antenna was modeled. The dipole is also larger than the monopole, since the smallest object to be modeled (the antenna) controls the uniform grid size in X and Y directions and size limitations of the FDTD code. The details about the grid size and the geometry of the soil medium for this case are tabulated in Table 3(a).
To decrease the computation cost, a much larger grid cell (3 cm in X and Y directions, and 5 cm in Z direction) was modeled (Table 3(a)). To satisfy sampling limitations (grid size < λ / 10) and study the scaling effect, the wavelength should be larger. Therefore, the frequency was selected to be 100 MHz (lower than 1.5 GHz in Case A). To satisfy the Courant’s condition for the new grid size, the time increment was increased to Δt = 50 psec.
|Simulated grid||149 × 149 × 69|
|Grid cell size||3 cm × 3 cm × 5 cm|
|Entire grid size||444 cm × 444 cm × 340 cm|
|Soil thickness||305 cm|
|Air thickness||35 cm|
The soil medium is exactly the same fully water-saturated sandy soil modeled in the previous case with 17% gravimetric moisture content. However, dielectric properties of the dispersive soil at the different frequency and time increments (f = 100 MHz, and Δt = 50 psec) are different. Therefore, the dielectric constant and coefficients (εAve, a1, b0, b1, and b2) of the Z-transform function required to model the dispersive electrical conductivity of the soil were recomputed for the new frequency and time increment. The new soil parameters are listed in Table 3(b). A center-fed resistively tapered ½ wavelength dipole antenna is modeled as the transmitter. The particular details of the resistive dipole are avoided by modeling the antenna electromagnetically as simply a tapered half-wave surface current source residing on the exposed coaxial insulator (maximum at the center, the point where the feed line joins the elements, and zero at the ends of the elements). This type of antenna may be used in a PVC-lined borehole filled with water. Therefore, the model simulates the antenna surrounded by water. Obviously, to model the dispersive nature of water and maintain the symmetry and accuracy on the circular interface around the antenna, water is modeled using the same technique used to model lossy dispersive soils (Weedon & Rappaport, 1997). For the same reason, the dielectric portion is modeled using the same technique used for lossy dispersive soils, despite the non-lossy and non-dispersive nature of the dielectric material. The PVC casing was ignored during the simulation to simplify the
Due to solving the problem at Δt = 50 psec, the FDTD code is very sensitive, and all digits are necessary to satisfy the stability conditions.
modeling and because the wall of the PVC is very thin compared to the wavelength (780 mm) of the EM wave. The dipole antenna is Z-polarized and the excitation signal is a 100 MHz cosine-modulated Gaussian pulse, progressively delayed along the antenna in the Z-direction (i.e., points along the Z-directed dipole are excited with a progressive phase delay proportionate to the traveling time of the current fed through the midpoint and along the dielectric portion of the dipole). Table 3(c) summarizes the details about the structure of the dipole antenna.
|Antenna depth||1800 mm|
|Borehole diameter||240 mm|
|Perfectly conducting core wire thickness||22 mm|
|Extended dielectric length||500 mm|
|Extended dielectric thickness||64 mm|
|Perfectly conducting outer conductor (shield) thickness||43 mm|
|Depth of water-filled borehole||500 mm|
|Gaussian width||10 nsec|
|Gaussian peak||75 nsec|
The dielectric constant and effective electrical conductivity of the extended dielectric are respectively assumed to be 2.1 and zero (Ω-1).
First, the wave propagation through the soil background was analyzed. Then, a rectangular DNAPL pool was modeled within the soil medium. Fig. 7 schematically shows the geometry of this DNAPL pool. The details about the geometry of the DNAPL pool scatterer are listed in Table 3(d).
|DNAPL Pool Geometry||Size|
|Horizontal area||45 cm × 45 cm × 15 cm|
|Clear distance to the antenna||22.5 cm|
|Coordinate of the pool center*||90 cm, 0 cm, 45 cm|
With respect to the center of the grid
Similar to Case A, the transmitting antenna is modeled in the code, but rather than modeling receiving antennae, the three different components of the electric and magnetic fields are computed at all grid points on the following depth and cross-sectional slices.
A cross-sectional (horizontal: XY-plane) slice, cutting through the antenna and DNAPL pool at the depth of 90 cm. Z and X components of the electric field (Ez and Ex) are shown on this slice (Fig. 8).
A depth (vertical: XZ-plane) slice, passing through the antenna and DNAPL pool. The magnitude of the power, derived from both Ex and Ez, is shown on this slice in Fig. 9 (Ey is zero on this slice due to symmetry).
As before, for each slice and each electric field component, three figures are shown: (i) incident (background) field, (ii) total field, and (iii) scattered field. As seen in Figs. 8(a) and 8(d), background results are symmetric. The total fields of Figs. 8(b) and 8(e) and the scattered field shown in Figs. 8(c) and 8(f) are asymmetric. The interesting and encouraging point is the visibility of the DNAPL pool over the entire medium within the total and scattered fields, even on the far side of the pool from the antenna. This predicts that dipole antenna boreholes of the CWR method can be drilled far from DNAPL-contaminated zones to reduce the risk of further vertically downward penetration of DNAPLs associated with drilling through contaminated zones. This appears to be more valid at higher degrees of water-saturation. Again, this potential can also be demonstrated in a different form as shown in Fig. 9.
Previously, in the case of the monopole transmitter, the received total and incident signals were computed at a receiver located right above the DNAPL pool. Now, the two are computed for a receiver located 175 cm above the pool to examine the possibility of minimizing the destructive effect of placing the receiving antenna too close to the pool. This figure again demonstrates a strong magnitude and travel time difference between the total and incident fields received at the receiver located far above the DNAPL pool. The strong magnitude difference (around 40%) and time difference (around 2.5 nsec.) between the two signals, once again, embraces the potential of the use of the cross-borehole GPR method to detect DNAPL pools. As in the monopole case, the early arrival of the total field and its higher magnitude can be respectively justified by the higher velocity of EM waves through the DNAPL pool due to its lower dielectric permittivity (relative to the water-saturated soil) and the lower loss due to the pool’s lower electrical conductivity. The dry soil case has been studied (it does not fit in the extent of this chapter) and proved to have weaker scattered signals, but still strong enough to have the potential to detect DNAPL pools using relatively widely spaced antennae.
As expected for the radiation pattern of dipoles, most of the energy is transmitted perpendicularly to the antenna through the mid-part of the antenna into the soil. The modeled dipole antenna is Z-polarized, and thus the Z-component of the electric field is the major component. The monopole antenna behaves similarly to the dipole antenna, the only difference being the strong signature of the DNAPL pool on the Z-component as well as X and Y components in the case of the dipole antennae. This wider spread perturbation due to the clutter promises more potential detection key-points for the dipole antenna installation. It is noteworthy that the perturbation due to the object on the X and Y components of the electric field of both monopole and dipole antennae appears more to the sides of the contaminated zone and perpendicular to the line connecting the center of the contaminated zone and the antennae. The perturbation on the Ez component in the monopole antenna case is distributed throughout the contaminated zone, while for the dipole one, this perturbation spreads to the far side across the contaminated zone as well.
6. Comparison with monopole experimentation
In this section, the numerically simulated monopole case is compared with experimentation for validation. The experimental setup (for more information refer to Farid et al., 2006) uses two PVC-cased ferrite-bead-jacketed monopole antennae connected to a vector network analyzer (Agilent 8714ES), and frequency-response measurements were collected for a homogeneous water-saturated sandy soil background. Fig. 11 shows a schematic of the experiment.
The frequency-response measurements are collected across the frequency range of 0.4 to 2.2 GHz. The frequency-response measurements are, then, transformed to the time domain using an inverse fast Fourier transform (IFFT) and an assumption of a narrow-width, wideband Gaussian pulse as the transmitted signal. Both the experiment and the FDTD model use the same Gaussian pulse source. Due to the frequency range used in the experimentation (0.4 GHz to 2.2 GHz), the width of the Gaussian signal should not exceed a maximum of: Wmin. = 1 / fmax. = 1 / 2.2 GHz = 0.455 nsec. Using narrower signals seems to create a noise. Therefore, the upper limit (0.455 nsec.) is selected as the width of the signal. Fig. 12 shows the transmitted signal.
The FDTD-simulated received signal at the bottom of the receiving antenna is shown in Fig. 13(a). As seen, this received signal is distorted and does not resemble the Gaussian transmitted one. Therefore, its peak is not easily distinguishable, since the received signal is modulated and noisy. To resolve this issue, the received signal should first be demodulated and then low-pass-filtered. The demodulation frequency can be found by observing the received signal in the frequency domain (computed via a fast Fourier transform). This demodulation frequency is observed to be dependent on the separation between the transmitting and receiving antennae. A MATLAB code was prepared to automatically observe the received signals in the frequency domain, find the proper demodulation frequencies, and find the proper low-pass filter to filter the noise. The processed received signal is shown in Fig. 13(b).
To transform the experimental frequency-response to the travel-time, the transmitted signal is fast Fourier transformed to the frequency domain and multiplied by the frequency-response (Fig. 14) to find the received signal at the receiver in the frequency domain. Then, the result is transformed back to the time domain using an inverse fast Fourier transform. The result (received signal in the time domain) is shown in Fig. 15(a).
This signal does not resemble the transmitted Gaussian signal. Therefore, it needs to be processed (demodulated and low-pass-filtered). The processed received signal is shown in Fig. 15(b). The demodulation frequency and filter design vary with the distance between the transmitting and receiving antennae, which is automatically calculated using the above-mentioned MATLAB code.
Since the vector network analyzer is calibrated at the end of the cables (the connection points to the monopole antennae), the experimental travel-time is measured between these two points. On the other hand, the FDTD travel-time (t3 –t1) is computed between the feed
cable and the bottom of the receiving antenna. Therefore, it needs to be adjusted for this difference.
Up to this point, the FDTD travel-time (t3 – t1) from the feed cable to the tip of the receiving antenna is computed. The travel time through the receiving antenna (t4 – t3), which is by symmetry equal to (t2 – t1), should be added to (t3 – t1) to find the total travel time between the feed and receiver cables (t4 – t1) for the FDTD model. The resulting travel time from the FDTD simulation can be used for comparison with the experimental results.
The travel time computed from the forward model is (4500 + 900 - 1000) × 2 psec = 8.8 nsec, which closely agrees with the one indirectly computed from the experimentally collected frequency-response data: (5700 - 1000) × 1.87 psec = 8.6 nsec. The difference is due to the slight, potential discrepancy between the dielectric constant assigned to the forward model (used from the results of another work by the authors (Zhan et al., 2007)) and the real values of the experimentation.
The intensities of the unprocessed received signals from the FDTD simulation (Fig. 13(a)) and experimentation (Fig. 15(a)) agree relatively well, but not perfectly. The reason is the potential slight discrepancy between the electrical conductivity assigned to the FDTD model compared to the actual one of the experiment. However, due to the difference between the necessary processing methods (different filters), the intensity of the processed received signals for the FDTD simulation (Fig. 13(b)) and the one of the experiment (Fig. 15(b)) do not agree as closely.
This comparison consists of the incident field for the homogeneous background soil. The comparison for the total and scattered fields at the presence of any anomalies (e.g., dielectric objects) will be conducted in the future.
A finite difference time domain (FDTD) model was developed for monopole and dipole antennae. Then, the scattering due to dielectric materials (to simulate DNAPL pools) in soils was modeled and analyzed. Results of the two simulated cases using the FDTD model demonstrate strong perturbation by the DNAPL pool on the electric field in the fully water-saturated sandy soil. In the case of the monopole antenna, the DNAPL pool target is more visible on the X and Y components of the electric field compared to the major component Z. The perturbation on the intensity of the electric field (|E|) transmitted by the monopole antenna is not as strongly visible as in the dipole case. In the dipole case, X and Y components are those parallel to likely hydraulic-conductivity contrast planes (e.g., usually horizontal clay lenses within a thick sand layer), which are potential locations to accumulate DNAPLs.
Different components of the electric field can selectively be collected using receiving antennae with different polarizations from the polarization of the transmitting antenna (e.g., a horizontally-polarized receiving monopole antenna and a vertically-polarized transmitting monopole antenna). Therefore, designing the receiving antenna alignment and polarization to selectively collect electric field components parallel to a possible DNAPL pool may help to compensate for a stronger perturbation on the minor components (X and Y) of the electric field emitted from a Z-polarized monopole antenna. These minor components should be of a high enough signal to noise ratio.
In the case of the dipole antenna, all three components of the electric field in the fully water-saturated soil have almost equal detection potential. In both of the above cases, there is a strong dielectric contrast between the DNAPL pool and the water-saturated soil. However, different radiation patterns of the dipole antenna compared to the monopole antenna may make the dipole antenna more desirable for DNAPL detection.
Field problems can be scaled down in size along with scaling up the frequency in non-dispersive soils to achieve the proper geometry and frequency for simulation purposes. This linear scaling of frequency and size may not work as well for dispersive soils, since frequency-dependent dielectric properties of dispersive soils add nonlinearity to the scaling problem. Other conclusions follow.
Images provided by such simulations show the field distribution that exists throughout the subsurface (i.e., similar to filling the entire volume with receiver antennae), but the field can only be observed practically by placing a reasonable number of receiving antennae at key underground positions with the appropriate polarization. This research can be used to find the radiation patterns of different antenna types and the interaction of the radiated field with soil heterogeneities, which leads to a better understanding of subsurface wave behavior at these key positions and aids the selection of optimum antenna patterns to cover these key positions.
While the depth of contamination is a problem for surface-reflection methods (e.g., GPR), there are no theoretical depth limitations for CWR, except practical drilling limitations and cost. The separation limitations between transmitting and receiving antennae used for CWR still exist. However, CWR has the advantage of using a one-way traveling path (transmission), unlike the two-way traveling path of surface-reflection GPR. In addition, the strong reflecting air-soil interface in the surface-reflection GPR technique is eliminated in the CWR technique and replaced with a better-controlled coupling between the borehole antennae and surrounding soil.
The perturbation due to the DNAPL target is stronger for the greater dielectric permittivity contrast between DNAPL pools and highly moist soil, as opposed to DNAPL plumes with low DNAPL saturation and dryer soils.
The signal to noise ratio of the scattered field by DNAPL pools should be high enough for measurements. As seen in the figures, the scattered field is comparable to the incident field. Therefore, if the signal to noise ratio of the incident field is high enough for measurement, the scattered field will probably have a large enough signal to noise ratio to be measurable as well.
The results of this forward model with monopole and dipole antennae show that the field perturbation (scattered = total - incident) for relatively large DNAPL pools at high enough DNAPL saturation, is of the same order of magnitude as the incident signal. This proves DNAPL detection using CWR in water-saturated soils feasible. The simulation tool can also be used as a forward model to develop an inverse scheme for DNAPL imaging.
Armed with the background data as well as the radiation patterns of different antennae (via simulations like those in this chapter), the existence of DNAPL pools can be confirmed with efficient inverse models and judicious placement of receiving antennae (i.e., pattern of antenna installation) where stronger perturbation and reception by receiving antennae are expected.
CWR may be a feasible and reasonable method to monitor DNAPL pools in a suitable environment. This most suitable environment is a medium consisting of a low-loss, low-heterogeneity porous material. In other media, it is more difficult to distinguish DNAPL accumulation from geologic variations, which are more complicated due to heterogeneity. Nevertheless, soil heterogeneity may not pose a crucial problem under water-saturated conditions since different soils behave similarly at relatively high degrees of water-saturation and high frequencies (the case is different for low frequencies). Monitoring DNAPL movement may well be possible or easier in an even less saturated heterogeneous environment because of the static nature of stratigraphic events and the dynamic nature of DNAPL flow. Several features of DNAPL pools may help to distinguish them from stratigraphic events, such as their irregular shapes with sharp lateral boundaries.
Finally, the FDTD model was compared for the incident field due to the monopole case in a homogeneous water-saturated sandy soil background with the experimental results. The reasonable agreement between both the travel time and intensity of the unprocessed, simulated and experimental results validates the FDTD model. The comparison and validation for the total and scattered fields at the presence of any anomalies (e.g., dielectric objects) need to be studied in the future.
This research was supported in part by the Gordon Center for Subsurface Sensing and Imaging Systems (CenSSIS), under the Engineering Research Centers Program of the National Science Foundation (NSF: Award Number EEC-9986821).
The authors would like to express gratitude for financial and scientific support provided by the Gordon CenSSIS and NSF.