Time-lapse electrical resistivity anomalies due to contaminant transport around landfills

The extent of landfill leachate can be delineated by geo-electrical imaging as a response to the varying electrical resistivity in the contaminated area. This research was based on a combination of hydrogeological numerical simulation followed by geophysical forward and inversion modeling performed to evaluate the migration of a contaminant plume from a landfill. As a first step, groundwater flow and contaminant transport was simulated using the finite elements numerical modeling software FEFLOW. The extent of the contaminant plume was acquired through a hydrogeological model depicting the distributions of leachate concentration in the system. Next, based on the empirical relationship between the concentration and electrical conductivity of the leachate in the porous media, the corresponding geo-electrical structure was derived from the hydrogeological model. Finally, forward and inversion computations of geo-electrical anomalies were performed using the finite difference numerical modeling software DCIP2D/DCIP3D. The image obtained by geophysical inversion of the electric data was expected to be consistent with the initial hydrogeological model, as described by the distribution of leachate concentration. Numerical case studies were conducted for various geological conditions, hydraulic parameters and electrode arrays, from which conclusions were drawn regarding the suitability of the methodology to assess simple to more complex geo-electrical models. Thus, optimal mapping and monitoring configurations were determined. Mailing address: Dr. Monica Radulescu, Conestoga Rovers and Associates, Waterloo, Ontario, N2V 1C2, Canada; e-mail: mradulescu@craworld.com


Introduction
Characterizing and monitoring the underground conditions of landfills and the location of subsurface contaminants has always been a challenge.Conventional environmental moni-toring used to determine the spread and fate of groundwater contamination is performed by the expensive and labour intensive task of drilling closely spaced boreholes for point sampling (Granato and Smith, 1999).More advanced, non-invasive geophysical techniques, less costly and more environmentally friendly, were developed as tools for interpreting the underground structures without disturbing them.Contaminant plumes usually have a sufficiently high contrast in physical properties against the host media due to an increase in dissolved salts in the groundwater and a resulting decrease in pore water resistivity; therefore, they may be detected by geo-electrical techniques.
Early studies have shown that the electrical resistivity method is one of the most efficient geophysical tools in detecting, delineating and monitoring groundwater contamination (Benson et al, 1983;Greenhouse and Slaine, 1983).More recently, Khair and Skokan (1998) investigated the early detection of salt water intrusion through systematic observation of electrical resistivity in selected positions with fixed electrode arrays.Geo-electrical techniques were also used in the detection of gasoline and saline leaks (Bevc and Morrison, 1991;Ramirez et al., 1996a), characterization of hazardous waste disposal sites (Ogilvy et al., 1999;Simmons et al., 2002), and the mapping and monitoring of contaminant plumes in groundwater systems (Ramirez, 1996b,c;Tekzan, 1999;Naudet et al., 2003).Electrical resistivity tomography was used to determine the spreading pattern of a saline tracer into groundwater (Singha and Gorelick, 2005).
However, field geo-electrical experiments are very expensive.In addition, characterizing landfills and location of subsurface contaminants is problematical as the extent and location of contaminant plumes are largely unknown during real site investigation.In order to minimize site characterization costs, it is advantageous to investigate theoretical responses of planned experiments before field programs are initiated.In addition, geophysical forward studies followed by inversion are essential to the final interpretation of observed geophysical data.
Electrical conductivity governing the electrical field is the key physical parameter used in geo-electrical exploration.As the dissolved contaminant travels with groundwater flow, advection and hydrodynamic dispersion cause the spreading of the contaminant plumes in directions transverse and longitudinal to the flow path.Consequently, the electrical conductivity of the targets varies not only with space but also with time.Thus it is essential to include the time-dependent feature of electrical conductivity in the forward and inversion studies performed during geo-electrical exploration.
The present research used the hydrogeological modeling of a groundwater contaminant flow path and the resulting distribution of resistivity/conductivity as a basis for the geophysical modeling of the geo-electrical structure.Images of the resistivity/conductivity distribution acquired through forward geophysical computations followed by inversion were expected to predict the initial distribution of concentrations in the hydrogeological model.The presented models and scenarios were designed to conceptualize a diverse range of hydrogeological environments, from simple to more complex, with various types of deposits and shapes of the plume.

Methodology
The methodology used in the paper has been developed by Yang (2005) for the theoretical investigation of geo-electrical responses associated with hydrothermal fluid circulation in a mid ocean ridge environment.The study presents the hypothetical, highly conceptualized scenario of a landfill, designed to assess the suitability of geo-electrical methods in describing the spreading of leachate plumes into the geologic environment.Several hypothetical scenarios of the leachate spreading into different geological structures were investigated.
The hydrogeological model simulated groundwater flow and contaminant transport by using FEFLOW 5.1, a 3D finite element code developed by WASY (WASY, 2005) capable of performing numerical modeling of density-dependent fluid flow and mass transport based on the hydrogeological components included in the conceptual model.
Once the spatial-temporal pattern of contamination was obtained by hydrogeological modeling as a distribution of concentrations, the concentration data was converted into resistivities/ conductivities, depicting the geo-electrical model.The conversion was based on the empirical Archie's law, describing the relation between concentration of contaminant and resistivity of groundwater.
The resulting electrical anomalies were computed by geophysical electrical modeling, using a finite difference software DCIP2D/ DCIP3D (2004), with different electrode configurations.
Electrical resistivity profiles and pseudosections corresponding to various electrode configurations were obtained by forward computation of the geo-electrical model.The direct conversion of salinity data acquired over a finely discretized mesh into electrical response offered the possibility of simulating the use of very large numbers of locations for the electrodes.Based on the results obtained at this stage, the configurations of electrodes which described most accurately the transient and space dependent image of the geo-electrical structure associated with the contaminant were selected.Geophysical inversion of the electrical data was then performed in order to recover the initial geo-electrical models.The inverted geo-electrical models were compared with the electrical models obtained previously by hydrogeological simulation.
Further analysis was conducted to evaluate the most suitable electrode configurations able to describe the spatial and temporal evolution of the contaminating plume originating from the landfill and consequently, the viability of the proposed methodology.

Hydrogeological modeling
Flow modeling is based on Darcy's law in porous media.As a density-dependent transport phenomena is simulated, FEFLOW uses the extended Oberbeck -Boussinesq approximation (Nield and Bejan, 1999).The solution of the governing equations is achieved using an implicit adaptive time stepping scheme.FEFLOW performs calculations in discrete time steps, imposed by the stability criteria required during the numerical simulation.Numerical stabilization is achieved with Petrov-Galerkin leastsquare upwinding, an alternative numerical scheme used to solve advective dominant flow and transport (Nguyen and Reynen, 1984).
A two-dimensional transient ground-water flow model was developed, describing various geological settings, according to different hydraulic conditions.The hydrogeological models are displayed in fig.1a-c.The hydrostratigraphical conceptualization of the models was chosen to depict gradually more complex geological and hydraulic settings in order to assess the viability of the proposed methodology on a wide range of conditions.
Triangle, a specialized code for creating twodimensional finite element meshes (Shewchuk, 1996(Shewchuk, , 2002)), was used to build Delaunay triangulations during mesh generation.A minimum angle of 20 degrees was used for each triangle, aiming to obtain a more uniform structure of the mesh while supporting the structure of the mod- el.A triangular finite elements mesh with dimensions of 600 m by 50 m, consisting of 5000 elements, was designed.The mesh was further refined in critical areas of contacts between different hydrostratigraphic units, to account for the numerical oscillations which could occur as a result of the increased hydraulic conductivity contrast.
The model simulates a transient evolution of the system over a period of 10000 days for all the presented models, steady state being usually reached after 2000-4000 days, depending on each model.The postprocessing of simulated data within FEFLOW was used as an evaluation tool for mass distribution and as a graphical output for fluid flow and mass transport.
The simulations were run under saturated conditions.The flow and the mass transport boundary conditions are presented in table I.The initial fluid velocity and contaminant concentra-tion throughout the model have been assumed to be equal to zero.A leachate concentration of 20400 mg/l was assigned over the interval BC on the upper boundary (fig.1a-c).The concentration of the leachate release was estimated as a conservative value based on the average composition of a waste landfill, consisting of specific chemical substances (McBean et al., 1994;Timur et al., 2000).The Total Dissolved Solid (TDS) values are converted to a common unit, NaCl and the unit expressed in terms of concentration (mg/l).At this stage of the study, focussing on finding the best electrical configurations which describe accurately the leachate spreading into the geological environment, the authors did not account for the chemistry of the plume, as it was not crucial for the purpose of the simulation.However, given the high concentration of the contaminant mass, simulations were run accounting for the density effect of the spreading plume.Klayer denotes the hydraulic conductivity of the intervening layers and Khost the hydraulic conductivity of the host rock.Porosities are regarded as φlayer, φhost and φfault.
Table II presents a summary of the hydraulic parameters used for the simulations.
Model 1 (fig.1a) presents the simple case of one intervening layer of lower hydraulic conductivity.
Model 3 describes more complex geological settings, with the presence of a layer displaced by a dip slip, normal fault, dipping in the groundwater flow direction (fig.1c).
Porosities and hydraulic conductivities have been assigned for the various conceptual scenarios according to the type of rocks and values taken from literature (Domenico and Schwartz, 1990).Less permeable layers were assigned lower hydraulic conductivities of 10 -6 m/s to 10 -8 m/s, characteristic to sandstones, while more porous layers were assigned higher conductivities of 10 -5 m/s, typical of gravel aquifers.
Over the simulated models and scenarios, representative contaminant plume distributions corresponding to specific time steps were selected to be converted into electrical data.The following time steps were selected for further geophysical forward and inversion computation: 2000 days for Model 1 (one intervening layer of lower hydraulic conductivity), 700 days for Model 2, Scenario 1 (discontinuous layer of higher hydraulic conductivity), 700 days for Model 2, Scenario 2 (discontinuous layer of lower hydraulic conductivity) and 300 days for Model 3, the faulted system.

Changes in bulk resistivity caused by contaminant plume infiltration
The mechanism generating the variations in resistivity is the change in fluid concentration created by the highly conducting contaminant fluid.Archie's law formula (4.1) was used to estimate the bulk resistivity of the rock as a function of the pore fluid's resistivity.where ρ f resistivity of the pore fluid, φ porosity, a and m empirical parameters set up for moderately well cemented sandstones, a=0.62 and cementing factor m=1.72.
According to various authors (Keller and Frischknecht, 1966;Levannier and Delhomme, 2003;Hwang et al., 2004), different formulas have been found in the literature, relating fluid concentration to bulk resistivity.The equation used in this study (4.2) (Morrison and Becker, 2004) has been derived from the resistivityconcentration plots (Keller and Frischknecht, 1966 where c= concentration of the fluid [mg/l]. The resistivity -concentration relationship was tested against various data sets obtained by direct resistivity measurements performed on a wide variety of saline solutes.Further verification showed that eq. ( 4.2) works the best at concentrations exceeding 3000 mg/l, which is the case for the landfill leachate in the present research.
The plume spreading pattern as concentration was converted directly into a model of resistivity/conductivity as a means of comparison with the resistivity model generated later during the geophysical forward/inversion procedure.According to the relationship between resistivity and concentration (eq.(4.2)), the distribution pattern of resistivities was similar to the concentration model.

Building up the input files for the DCIPF 2D software
The concentration data corresponding to the nodal values of the hydrogeological model built on a finite triangular elements mesh was discretized into rectangular distribution and, consequently, converted to time-dependent resistivity/conductivity models.The F2D converter software was used to build up the input files for the forward geophysical modeling based on the electrical resistivity data.Various arrays with different electrode configurations were used to derive potentials from the resistivity data, as follows: pole-pole (pp), pole-dipole array with the potential electrodes on the left (pdL), pole-dipole array with the potential electrodes on the right (pdR), dipole-dipole (dd).
The forward modeling of the direct current DC potentials computed time-lapse apparent resistivity anomalies based on a finite difference technique used by the DCIPF2D code (DCIP2D, 2004).The anomalies consisted of synthetic data that would be acquired over the 2D resistivity structures obtained previously from the concentrations acquired through hydrogeological modeling.
Initially, the two dimensional, vertical cross sectional domain was set 600 m wide and 50 m deep, corresponding to 60 cells on the horizontal dimension and 10 cells on the vertical.The model was extended according to the mesh generation algorithm (Li et al., 1995) to 620 cells in the x-direction (x=−2042 ... +2642 m) and 54 cells in the z-direction (z=0 ... −170 m).For each configuration, surface electrodes were spaced gradually in incremental steps of 1 to 5 m in the interval x=(0 ... +600) m, until reaching the right side of the model.The initial data was contaminated with independent Gaus-sian noise whose standard deviation was equal to 5% of each accurate datum.
The resulting observations of forward modeling consisted of images of apparent resistivity distribution as pseudosections acquired with different configurations of electrodes over the initial geo-electrical model.
Additional information was provided by plots of the apparent resistivity variation over the surface of the model.

Inversion of the DC data
The DC potentials as apparent resistivities were inverted to recover the geo-electrical model by using the computing program DCINV2D (DCINV2D, 2004;Oldenburg et al., 1993).
The same mesh from the forward simulation was kept for the inverse problem, such that the inversion generated resistivities for all the 620 horizontal cells until the observations were adequately fit.
The inversion program estimates a reference model and an initial model, assigns estimated error standard deviations to the data, and constructs a smooth model that best fits the data by reproducing the initial electrical model to the expected value.

Results and discussions
The figures show the results of the research by describing the evolution of the contaminant plume for various geological models, as follows: the evolution in space and time of the solute transport, in terms of concentration, as acquired during the hydrogeological numerical modeling generated with FEFLOW (figs. 2 to 5); the associated geo-electrical response, based on the geophysical forward modeling performed with the DCIPF2D code.The geo-electrical response is described in terms of pseudosections (figs.7a-d to 10a-d) and apparent resistivity profiles (only for Model 2, Scenario 1, fig. 6).Figures 11a-c to 14a-c present the result of the geophysical inversion through the DCINV2D code, as a summary of observed data, predicted data and difference between observed and predicted data normalized by standard deviation.Figures 15 to 18 show the recovered resistivity models.
Finally, the computed resistivity models were compared with the models describing the hydrogeological mass distribution.

Hydrogeological modeling
For Model 1, the plume was constrained by a layer with a hydraulic conductivity two orders of magnitude lower than the host rock.As a result, the contaminant is mainly contained above the lower hydraulic conductivity layer (fig.2).
The presence of intervening layers and a fault, with different hydraulic parameters and an increased contrast of hydraulic conductivities between different formations enables a more com-     plex pattern of the plume (Freeze and Witherspoon, 1967).As a result, the next conceptual models relied on a higher contrast of hydraulic properties between intervening layers, host rock and fault.A higher conductivity of the host rock enables a more diffuse spreading of the contaminant, at lower concentrations (fig.4) compared to a lower hydraulic conductivity of the host rock, which favours the development of a more sharply delineated plume, at higher concentrations (fig.3).
For Model 2, Scenario 1, discontinuous layer, two orders of magnitude contrast between the higher conductivity layer and the medium (K layer = 0.125 * 10 −4 m/s, Khost= 0.00042 * 10 −4 m/s) can be regarded as a physical shield, causing refraction of the flow line such that flow in the higher conductivity layer is mainly horizontal, meanwhile flow in the lower conductivity medium in essentially vertical (Freeze andWitherspoon, 1967, Neuman andWhiterspoon, 1969).As a result of this phenomenon, the contaminant plume spreads at higher concentrations mainly above and below the higher conductivity layer (fig.3).
Figure 4, Model 2, Scenario 2, describes the same discontinuous layer model, except its hydraulic conductivity is two orders of magnitude lower than the hydraulic conductivity of the host rock (K layer = 0.00042 * 10 −4 m/s, Khost = 0.125* *10 −4 m/s).The layer acts as a barrier to the plume development.Significant flow occurs mainly within the higher conductivity rock, enabling diffuse transport of the plume in the system, at lower concentrations (fig.4).
A comparison between Scenario 1 and Scenario 2 shows that the presence of a lower conductivity layer (Scenario 2) (fig.4) favoured the retention of high concentrations of contaminant in a small area, meanwhile the overall spreading of the contaminant in the host rock at lower concentrations was significantly enhanced compared to the case of the higher hydraulic conductivity layers (Scenario 1) where the contaminant spread as a more sharply delineated, elongated plume of higher concentrations (fig.3).
Model 3 describes a dip slip, normal fault, dipping in the same direction as the fluid flow.The hydraulic conductivities were 1.36 * 10 −4 m/s  Overall, the density effect due to high concentration of the pollutant generated a decrease in fluid flow velocities.A higher hydraulic conductivity layer associated with reduced fluid velocity enabled the accumulation of the contaminant within the layer.The phenomena can be regarded as a vertical restriction of the plume development.

Forward geophysical modeling
For each array, the electrodes spacing was gradually increased incrementally until reaching the right limit of the model on the surface.The incremental spacing was varied from steps of 1 to 10 m at a time.As a preliminary observation, a lower incremental spacing between electrodes along the entire length of the model enables higher resolution of the apparent resistivity pseudosections.A major outcome of the proposed methodology resulted from the versatility of us-  ing a wide range of steps for electrode spacing, where field and technical constraints can be disregarded and consequently, detailed information about the subsurface can be depicted.It is a crucial aspect, as the proposed methodology is a theoretical exercise intended to offer the best preparatory knowledge before the investigative work in the field.In this case, the only constraints referred to the computing procedures, where significant running times restricted the use of very small incremental steps.The electrodes spacing which most accurately describes the horizon where the plume is moving within a reasonable computing time was selected.Without previous hydrogeological information, the decision of which electrode spacing is most appropriate is quite difficult to reach.
The electrical resistivity figures in the next sections present the case of an increment of 5 m in electrodes spacing.
Based on the concentration values transformed into resistivities, the geophysical for-Fig.13a-c.Result of geophysical inversion, pdL configuration, Model 2, Scenario 2, at 700 days: a) observed data; b) predicted data; c) difference between observed and predicted data.

Apparent resistivity profiling
The following observations were drawn based on the apparent resistivity profiling using various electrode configurations for each of the    hydrogeological scenarios.Figure 6 presents the apparent resistivity profiling for Model 2, Scenario 1, at 700 days.Apparent resistivities profiling cannot offer relevant information about the distribution of the resistivities with depth.However, the profiles can be regarded as useful to delineate the contact between the noncontaminated/contaminated area at the surface, along the investigation line.
The profiles showed decreased values of resistivity in the leaching area, increasing with time and space according to the spreading and dilution of the plume, until eventually reaching again the base line of the noncontaminated area.The initial flat trend of the apparent resistivity, corresponding to the presence of the noncontaminated area at designated intervals of time, can be used as a baseline for characterizing clean conditions.
However, there were certain resolution differences which made some configurations more suitable to delineate the leachate spread at the surface, as follows: -In the case of dd configuration, the initial increasing peak on the profile distinctively marked the contrasting properties between the noncontaminated and contaminated area, offering accurate information about where the leachate spill originated at the surface.Figure 6-a, for example, showed a sharp delineation of the initiation of the plume.According to it, the length of the noncontaminated area at the surface is of approximately 170 m, in fair agreement with fig.3-c, snapshot of concentrations distribution over the model.
-For the pdL configuration, the sharply increasing peak on the profile before the clean baseline, as a result of the contrasting properties between the noncontaminated/contaminated surface, indicated the end of the contaminated area and could be used as marker as long as the contamination does not extend outside the horizontal boundaries of the initial geo-electric model, which is the case in fig.6-b.
-For the pdR configuration, the initial increasing peak on the profile marked sharply the beginning of the contamination area, providing reliable information about the leachate spreading at the surface (fig.6-c).
-In the case of pp configuration, the peak present for other electrodes configurations was absent at both ends of the contaminate area (fig.6-d).
As a first conclusion, the apparent resistivity profiles can be used as a rougher interpretation tool to delineate the contaminant plume spreading through time at the surface.dd and pdR are the optimal configurations which delineate the initiation of the leachate spill at the surface.

Apparent resistivity pseudosections
Pseudosections of apparent resistivities generated by the leachate spreading were investigated with various electrode configurations (dd, pdL, pdR, pp), respectively figs. 7a-d, 8a-d, 9a-d and 10a-d.For Model 1, of a uniform environment and consequently sharp delineation of plume, all electrode configurations described quite accurately the extension of the contaminant in time and space (fig.7a-d).
In the case of Model 2, Scenario 1 (fig.8ad), where the plume becomes more diffuse at lower depths, accurate information was obtained for the dd and pdL configurations.
For Model 2, Scenario 2 (fig.9a-d) and Model 3 (fig.10a-d) and a diffuse pattern of the contaminant, relatively accurate information was provided by the pdL configuration.
Dispersion along the flow lines in the porous media creates dilution of concentration, where the resistivity contrast between the leachate and the clean groundwater decreases and consequently the geo-electrical image of the model is less accurate.As first order estimates, not all the dipole electrode configurations were equally suitable to offer reliable information about the spreading of the contaminant.
All the electrical imaging pseudosections depicted accurately the patterns of contaminant for sharply delineated plumes, with concentrations ranging between 10000 to 12000 mg/l, regardless of the electrodes configuration.
Where the plume spreads diffusely into the groundwater, at lower gradients of concentration, between 2000-4000 mg/l, the resistivities contrasts have not been sufficiently sharp to create an accurate forward electrical image (figs. 7a-d and 8a-d).The pseudosections which have best revealed the time-space transient feature of the contaminant plume given complex geological and hydraulic settings were acquired with the pdL configuration.
Prior to inversion, results plotted as conventional pseudosections enabled preliminary interpretations, as first order estimates of data quality, anticipating the spreading of the plume in terms of structural features of the geo-electrical model.
Based on the ability of the pdL configuration to profile the extent of the plume at the surface and to offer an acceptable pseudo image of the distribution of apparent resistivities in depth, the pdL was chosen to perform the inverse problem of the geo-electrical model.

Inversion of the DC data
The inversion was carried out for the selected configuration by supplying to the program the DC data observation file from the forward modeling and the same user-designed mesh.The inversion program estimated a reference model and an initial model, assigned estimated error standard deviations to the data, and constructed a smooth model that best fits the observed data to the expected value.Figures 11ac to 14a-c present for each inverted model the initial resistivity model, the observed data from the forward modeling, the predicted data, and the difference between the observed-predicted data normalized by standard deviation.As can be observed from figs.15 to 18, the essential features of the resistivity were all recovered.

Conclusions and further investigations
The hydrogeological numerical simulation of contaminants spreading in the groundwater was used as a support for the geophysical modeling of the associated geo-electrical models.Four electrode configurations, pole-pole array (pp), pole-dipole with the potential electrodes on the left (pdL), pole-dipole with the potential electrodes on the right (pdR), dipole-dipole (dd), were used for the numerical modeling of the electrical anomalies.The results were initially computed through forward modeling as apparent resistivity profiles and pseudosections, followed by inversion.
Reliable and accurate data from 1D profiles of electrical resistivity, delineating the initiation of the contaminating plume at the surface along the investigation line, can be acquired regardless the complexity of the model by the use of either dd or pdR configurations.The pdL apparent resistivity profiles have proven suitable to delineate the end of the contaminated area at the surface as long as the contamination did not extend beyond the limits of the investigation line.
The pseudosections, as preliminary interpretations prior to inversion, have revealed accurately the time-space transient feature of the contaminant plume for all the electrodes configurations, as long as the hydrogeological environment is fairly uniform and there is a sharp contrast of concentrations between the leachate and the fresh groundwater.For complex and more diffuse distributions of concentrations, with lower gradients between the resistivity features in the models, the pseudosections produced the best results for the pdL configuration.
However, the optimal contrast between resistivities cannot be quantified and generalized from one case study to another.The level of contamination detected by forward modeling of resistivity in different aquifers can be different even when using the same electrodes configurations and spacing, due to diverse geo-electrical structures and resistivity contrasts.
The inversion of the observed data recovered the essential resistivity features for all the scenarios.Consequently, there was a good time-space correlation between the recovered resistivities profiles and the plume as described by the hydrogeological model, for all the electrode configurations in case of the uniform model and for pdL configuration for more complex distribution of the plume.The inverted electrical structures also offered a matching distribution of the contaminant concentrations, as initially described in the hydrogeological model.
As reemphasized by the present research, a good knowledge of the contamination patterns in the area before planning field work in real situation is an essential asset.Based on the complexity of the geological and hydraulic settings, an educated decision can be made regarding the further use of certain array configurations, aiming to acquire the best results and optimization of costs.Within the selected array, the proposed methodology allows the theoretical use of a wide range of electrodes spacings, directly related to depth of investigation, from which the configuration offering the best information about the contamination horizon can be selected.
Furthermore, there is no limitation to the scale of the model, as runs can be performed for any dimensions of the finite elements and, therefore, finite difference meshes.

Fig. 2 .
Fig. 2. Evolution in space and time of the solute transport, in terms of concentration, Model 1.

Fig. 3 .
Fig. 3. Evolution in space and time of the solute transport, in terms of concentration, Model 2, Scenario 1.

Fig. 4 .
Fig. 4. Evolution in space and time of the solute transport, in terms of concentration, Model 2, Scenario 2.

Fig. 5 .
Fig. 5. Evolution in space and time of the solute transport, in terms of concentration, Model 3.

Fig
Fig. 11a-c.Result of geophysical inversion, pdL configuration, Model 1, at 2000 days: a) observed data; b) predicted data; c) difference between observed and predicted data.

Fig
Fig.12a-c.Result of geophysical inversion, pdL configuration, Model 2, Scenario 1, at 700 days: a) observed data, b) predicted data, c) difference between observed and predicted data.

Fig
Fig.14a-c.Result of geophysical inversion, pdL configuration, Model 3, at 300 days: a) observed data; b) predicted data; c) difference between observed and predicted data.

Table I .
Flow and the mass transport boundary conditions.

Table II .
Hydraulic parameters.