Flow dynamics in hyper-saline aquifers: hydro-geophysical monitoring and modeling

Saline–freshwater interaction in porous media is a phenomenon of practical interest particularly for the management of water resources in arid and semi-arid environments, where precious freshwater resources are threatened by seawater intrusion and where storage of freshwater in saline aquifers can be a viable option. Saline–freshwater interactions are controlled by physico-chemical processes that need to be accurately modeled. This in turn requires monitoring of these systems, a non-trivial task for which spatially extensive, high-resolution non-invasive techniques can provide key information. In this paper we present the field monitoring and numerical modeling components of an approach aimed at understanding complex saline–freshwater systems. The approach is applied to a freshwater injection experiment carried out in a hyper-saline aquifer near Cagliari (Sardinia, Italy). The experiment was monitored using time-lapse cross-hole electrical resistivity tomography (ERT). To investigate the flow dynamics, coupled numerical flow and transport modeling of the experiment was carried out using an advanced three-dimensional (3-D) density-driven flow-transport simulator. The simulation results were used to produce synthetic ERT inversion results to be compared against real field ERT results. This exercise demonstrates that the evolution of the freshwater bulb is strongly influenced by the system’s (even mild) hydraulic heterogeneities. The example also highlights how the joint use of ERT imaging and gravity-dependent flow and transport modeling give fundamental information for this type of study.


Introduction
Multiphase flow in porous media has been the subject of intensive study for many decades, motivated, amongst other factors, by important economic considerations linked to the petroleum industry. Another field where interaction of pore fluids having different physical properties, which is of particular importance, is saline-freshwater systems. In this case, important density and viscosity differences between saline and fresh waters control the relative motion and mixing of the two phases. Characterizing and modeling these coupled flow and transport phenomena is a very challenging task, particu-larly in the presence of the hydraulic heterogeneities always present in natural porous media (e.g., Werner et al., 2013;Ketabchi et al., 2016).
The most common situation where saline-freshwater systems have practical environmental and socio-economic implications is related to seawater intrusion in coastal aquifers, often exacerbated by overexploitation of groundwater, particularly in arid and semi-arid regions such as those surrounding the Mediterranean basin (e.g., Kallioras et al., 2010;Rey et al., 2013;Dentoni et al., 2015). Another context where the study of saline-freshwater interactions is highly important is the injection and storage of freshwater in brackish or salty aquifers for later use in agriculture or for domestic purposes, also known as aquifer storage and recovery (ASR; e.g., Pyne, 1995;Dillon, 2005).
Many studies of density-dependent flow and transport phenomena in porous media have been conducted over the past decades (e.g., Gambolati et al., 1999;Simmons et al., 2001;Diersch and Kolditz, 2002). Instabilities and fingering can take place when denser water overlies lighter water (e.g., Simmons et al., 2001). Ward et al. (2007) gave an introductive literature review on density-dependent modeling, with a particular focus on ASR. The first studies on the injection of freshwater into a saline aquifer were performed by Bear and Jacobs (1965) and Esmail and Kimbler (1967). The latter investigated the tilting of the saltwater-freshwater interface, a phenomenon known as "buoyancy stratification". More recent studies have analyzed the efficiency of ASR for both field and synthetic cases (e.g., Kumar and Kimbler, 1970;Moulder, 1970;Kimbler et al., 1975;Ward et al., 2007Ward et al., , 2008Lu et al., 2011;Zuurbier et al., 2014). Ward et al. (2008) conducted a numerical study to evaluate the efficiency of ASR under density-dependent conditions with anisotropy and heterogeneity of high and low permeable layers. Van Ginkel et al. (2014) studied the possibility of extracting saltwater below the freshwater injection to prevent the spreading of freshwater at the top of the aquifer. Alaghmand et al. (2015) investigated fresh river water injection into a saline floodplain aquifer and developed a numerical model for the optimization of injection scenarios.
The behavior of saline-freshwater systems becomes increasingly complex with larger density and viscosity contrasts. To date, very little research has been done on the effects of freshwater injection in highly saline aquifers that can reach total dissolved solids concentrations of 100 g L −1 . Understanding these complex systems is limited not only by the need to develop non-trivial coupled flow and transport models but also by the scarce availability of effective monitoring techniques. The latter are, under field conditions, typically limited to borehole measurements that can only provide point information in spatially heterogeneous hydraulic systems with time-changing salt concentrations.
As in many other subsurface characterization problems, a major contribution can be made by non-invasive, spatially extensive, geophysical techniques. In particular, electrical and electromagnetic methods are very suitable in the context of saline-freshwater interactions, since electrical conductivity varies over orders of magnitude depending on solute concentrations. While the use of these methods is common in seawater intrusion studies (e.g., Goldman and Kafri, 2006;Nguyen et al., 2009), only few studies have used geophysics to monitor ASR experiments. Davis et al. (2008) used time-lapse microgravity surveys to monitor the utilization of an abandoned coal mine as an artificial ASR site. Maliva et al. (2009) investigated the use of geophysical borehole logging tools applied to managed aquifer recharge systems, including ASR, to improve the characterization of aquifer properties. Minsley et al. (2011) developed an integrated hydro-geophysical methodology for the siting, operation and monitoring of ASR systems using electrical resistivity, time-domain electromagnetics and seismic methods. Parsekian et al. (2014) applied geoelectrical imaging of the subsurface below an aquifer recharge and recovery site alongside with hydrochemical measurements to identify preferential flow paths.
A major step forward in saline-freshwater systems monitoring can be made by improving the efficiency of advanced geophysical techniques, and electrical tomographic methods in particular. Electrical resistivity tomography (ERT) is widely used today in hydro-geological and environmental investigations. Often applied in tracer studies (e.g., Kemna et al., 2002;Vanderborght et al., 2005;Cassiani et al., 2006;Doetsch et al., 2012), ERT is a natural choice for salinefreshwater interaction monitoring, given the correlation between the salinity of a pore fluid and its electrical conductivity. Time-lapse ERT, where only the changes in electrical conductivity over time are imaged (e.g., Kemna et al., 2002;Singha and Gorelick, 2005;Perri et al., 2012), can be especially effective in tracking dynamic processes. Whereas tracer studies are typically designed with injection of a saline tracer into fresh surrounding groundwater, only very few studies have dealt with the inverse case of freshwater injection into a saline formation. For instance, Müller et al. (2010) conducted tracer tests also using a less dense tracer with lower electrical conductivity than the ambient groundwater, monitored with ERT.
The goal of this study is to present a general approach for the characterization, monitoring and modeling of complex saline-freshwater systems, based on the combination of noninvasive techniques and accurate numerical modeling. To our knowledge, no such a comprehensive hydro-geophysical approach concerning freshwater injection in saline aquifers has been presented so far in the scientific literature; thus, we believe this case study can be very useful as a starting point for other, more comprehensive methodological testing. In this study we limit ourselves to integrating field data and modeling in a loose manner, with no aim at this stage to develop a full data assimilation framework, as implemented elsewhere for simpler systems (e.g., Manoli et al., 2015;Rossi et al., 2015). The key message that can be derived from the joint use of advanced field techniques and advanced numerical modeling is nonetheless apparent in the presented case study, and more complete assimilation approaches are possible provided that the advantages and limitations of the individual components (data and models) are fully understood as shown in the present paper.
The approach is presented in the context of a case study where we injected freshwater into a hyper-saline aquifer in the Molentargius Saline Regional Park in southern Sardinia, Italy. The experiment was monitored using cross-hole timelapse ERT. To investigate the mixing processes, the resulting ERT images are compared with the results of a synthetic numerical study of the same experiment. We consider here both homogeneous and heterogeneous (layered) systems. For a quantitative comparison between the field and synthetic studies, spatial moments of the freshwater bulb are calculated.
2 Field experiment

Site description
The Molentargius Saline Regional Nature Park is located east of Cagliari in southern Sardinia, Italy (Fig. 1). The park is a wetland situated very close to the coastline. The exceptional nature of the site is given by the presence of both freshwater and salty water basins separated by a flat area with mainly dry features (called "Is Arenas"). The freshwater areas include two ponds that originated as meteoric water retention basins. The salty water areas include the stretches of water of the former system of the Cagliari salt pans.
The park area is characterized by an oligocenic-miocenic sedimentary succession of ca. 100 m (Ulzega and Hearty, 1986) overlaid by pleistocenic deposits of marine and continental origin and by alluvial and offshore bar deposits, whose origin is still debated Thiel et al., 2010). This ongoing scientific debate has implications for the comprehension of the phenomenon of hyper-saltiness of the park groundwater.
The specific site of investigation is located in the flat dry area within the park (Is Arenas, Fig. 1c). The water table of the unconfined aquifer is stable at 5.2 m below ground surface (b.g.s.), and practically neither lateral groundwater flow nor tidal effects are evident. The sediments are composed mostly of sands, with thin layers of silty sand, clayey sand and silty clay (Fig. 2). The groundwater reaches salinity levels as high as 3 times the NaCl concentration of seawater. Such high salt concentration is likely the long-term legacy of infiltration of hyper-saline solutions from the salt pans dating back, in this area, to Roman times. Electrical conductivity fluid logs (see Fig. 3) recorded in boreholes allowed two zones to be discriminated, with a transitional layer in between; (1) from the water table to a depth of 6.5 m the water electrical conductivity is about 2 S m −1 , and (2) below 12 m depth the water electrical conductivity reaches 18.5 S m −1 . Note that Fig. 3 also reports the time-lapse evolution of the vertical electrical resistivity profile as a result of the freshwater injection described in the following section.

Freshwater injection
Five boreholes for ERT measurements were drilled with 101 mm inner diameter to a depth of 20 m and positioned in the shape of a square with 8 m sides (four corner boreholes) and one borehole at the center (Fig. 1b). All boreholes are equipped with a fully screened PVC pipe (screen with 0.8 mm size).
In November 2011, 19.4 m 3 of freshwater with an electrical conductivity of 0.03 S m −1 , stored in a tank, was injected into the saline aquifer. This was done through the central borehole using a double packer system with an injection segment of 1 m length. The injection chamber was set between 13 m and 14 m b.g.s. The injection rate was entirely controlled by the natural pressure gradient, given by the water head in the tank and the depth of injection (i.e., 13 m to 14 m b.g.s. plus 2 m head in the tank above the surface). The natural pressure gradient provided for an initial injection rate of 0.5 L s −1 . However, during injection (after about 1.5 h) this rate immediately rose to a rate of about 2.75 L s −1 . We assume that this was due to a clogging of the backfill material, which was "de-clogged" after 1.5 h. In total, discharging the tank took about 4 h.

ERT monitoring
The direct electrical conductivity measurements described in the previous subsection correspond to the data that would be available as a result of a standard monitoring plan, and is highly insufficient for drawing any conclusions concerning the processes that take place during and after freshwater injection. The available dataset was great enriched by ERT measurements, described below.

Data acquisition
Time-lapse ERT monitoring was applied during the injection experiment in order to image the developing freshwater bulb, "visible" thanks to its lower electrical conductivity compared to the surrounding saltwater. Each borehole bears externally to the casing 24 stainless steel cylindrical electrodes, permanently installed from 0.6 to 19 m depth with 0.8 m separation, with the exception of the central borehole where the first electrode is placed at the surface and the last at 18.4 m depth. ERT measurements were carried out in a two-dimensional (2-D) fashion, along two vertical planes diagonal along the boreholes, i.e., one plane was using the borehole numbers 1, 5, and 3 and the second plane the borehole numbers 2, 5, and 4 (see Fig. 1b), thus making use of 72 electrodes per plane. This choice, in contrast to a full 3-D acquisition, was predicated on minimizing the acquisition time, given that the freshwater-saltwater movement was expected to be relatively rapid.
The ERT measurements were conducted using a Syscal Pro and adopting different configuration setups, consisting of in-hole dipole-dipole measurements in a skip-zero mode (i.e., adjacent electrodes form a dipole) and cross-hole dipole-dipole (hereafter referred to as bipole-bipole) measurements (Fig. 4). Measurements were collected in normal and reciprocal configurations (i.e., exchanging the current and potential dipoles) for estimation of data errors. The acquisition for one complete measurement frame (consisting of roughly 7300 individual readings) required about 40 min. . Schematic description of the ERT measurement configurations used. For dipole-dipole measurements, one dipole is always within one borehole, the other dipole also moves into the adjacent borehole. Bipole-bipole measurements are done as cross-hole measurements and are also changing as diagonals (i.e., A stays while B moves downwards for up to five electrode positions before A is also moved, similarly for M and N).
ERT data were acquired in a time-lapse manner to investigate the changes over time caused by the electrical conductivity changes of the developing freshwater bulb within the saline aquifer. The first time step, T 0, was acquired before the start of injection in order to compare the following individual time steps with the background image. These were measured on the day of injection, 1 day after injection, and 5 days after injection.

Data processing and time-lapse ERT inversion
Due to technical errors (such as bad connection of electrodes, problems with power supply) and varying data quality, the ERT data were processed prior to inversion. In particular, data having a misfit larger than 5 % between normal and reciprocal readings were removed.
The temperature difference between the groundwater (21 • C) and the injected freshwater (18 • C) was relatively small. Changes in electrical conductivity due to temperature effects are in this case about 5 % (see, e.g., Sen and Goode, 1992). Compared to the variation in electrical conductivity between the two fluids, which is about 3 orders of magnitude, the temperature effect is considered negligible.
The ERT field data from the freshwater injection experiment were inverted using the smoothness-constraint inversion code CRTomo. A full description of the code is given by Kemna (2000). In the inversion, the data errors are represented according to a linear model expressed as ε = a/R +b, where R is the measured electrical resistance. For the case at hand, the error parameters a (absolute) and b (relative) were set to 0.0001 m and 10 %, respectively.
Resistivity images exhibit a variable spatial resolution (e.g., Ramirez et al., 1995;Alumbaugh and Newman, 2000;Nguyen et al., 2009). A useful indicator for this variation is the cumulative sensitivity s (e.g., Kemna et al., 2002;Nguyen et al., 2009). The sensitivity indicates how a change in electrical resistivity of a certain model cell affects a transfer resistance measurement. Analogously, the cumulative sensitivity quantifies the change of a complete dataset to a changing model cell, and its analysis is an important step in the inversion process. Note that an objective choice for a threshold, which identifies zones where "reliable" vs. "unreliable" ERT imaging, is not feasible. In a more qualitative manner one can assume, empirically, that a cumulated sensitivity clearly below 1e − 3 leads to weak imaging. Figure 5 shows exemplarily the cumulative sensitivity distribution for the inversion of one dataset (image plane boreholes 1 − 5 − 3 at time T 0, i.e., the background image). The geometry of the boreholes and the electrodes, in combination with the employed measurement configurations, yields a relatively good coverage within the area of interest (i.e., mainly the area around the central borehole).
In a time-lapse monitoring framework, one is primarily interested in the temporal changes of data and parameters. Therefore, we used the "difference inversion" approach of time-lapse ERT (e.g., LaBrecque and Yang, 2000; Kemna et al., 2002), where the inversion results are changes with respect to the background data at time T 0. The advantage of this approach is that modeling errors and data errors correlated over time are canceled out to a significant degree and associated imaging artifacts that would occur in a standard inversion are suppressed.
K. Haaken et al.: Flow dynamics in hyper-saline aquifers

ERT imaging results
The ERT dataset was collected under challenging conditions, in particular as the very large salinity contrasts are manifested as extreme electrical conductivity differences over space and time. Large electrical conductivity can occasionally bring DC electrical currents into a nonlinear (non-Ohmic) regime, which in turn can lead to violation of the conditions for the reciprocity theorem (Binley et al., 1995;Cassiani et al., 2006). This has clear implications in terms of data processing, as in particular the error analysis based on reciprocal resistances may not guarantee that direct and reciprocal resistances are equal to each other. Filtering the data according to a reciprocity discrepancy equal to the data error level chosen for the inversion (see above) meant that a fairly large percentage of the data (about 50 %) were rejected. Nonetheless, a large volume of resistance data were still retained (nearly 2000 values per time instant).
The very high electrical conductivity of the system, which is characteristic of this experiment, has also another consequence; i.e., separated inversion of the different electrode configurations (dipole-dipole and bipole-bipole) showed that the bipole-bipole configurations provide better overall results than the dipole-dipole configuration results (not shown here). This is not a common situation, as observed elsewhere in situations of standard resistivity ranges (e.g., Deiana et al., 2007Deiana et al., , 2008, where dipole-dipole data provide higher-resolution images than bipole-bipole data that generally only give smoother images as information is averaged over large volumes. In the case shown here, for an in-hole current dipole, the current lines will not penetrate far away from the borehole as they are short-circuited by the large electrical conductivity of saline water surrounding at all times the external boreholes, whereas for the cross-hole current bipole the current lines "have to" penetrate through the volume between the boreholes. Thus, the sensitivity for the dipole-dipole configurations decreases very strongly with increasing distance from the boreholes. However, the dipoledipole configuration still manages to provide high sensitivity in the area close to the central borehole, particularly at measurement times where the freshwater bulb surrounds this borehole. Hence, the data coming from both configurations were used for inversion. Figure 6 shows the background image (time T 0) before the start of freshwater injection. The electrical resistivity of the saturated zone is very low and vertical changes due to layering of lithologies are not visible. Only a gradual change to higher resistivities in the upper part just below the water table can be seen. This can partly be attributed to the smoothness constraint applied in ERT inversion. However, this feature is also consistent with background conductivity logs (Fig. 3).
The obtained time-lapse ERT images of the freshwater injection experiment are shown in Fig. 7; the distribution of the injected freshwater in the aquifer surrounding the central borehole is clearly visible, in agreement with the time- lapse conductivity logs in Fig. 3. The very fast vertical migration of the freshwater plume is also apparent. Between 2 and 6 h after the start of injection, the injection borehole (and its surroundings) is nearly totally filled with freshwater, as confirmed by Fig. 3 (after 5 h). However, from the ERT images the freshwater also seems to move downwards below the injection chamber. A few hours after injection, the freshwater plume nearly disappeared in the ERT images, and 1 day after injection the ERT image seems to have gone back to the background situation (as also confirmed by the conductivity logs in Fig. 3).
At about 10 to 11 m depth, the difference images show a separation of the plume into two parts. A layer of finer sediments (see Fig. 2) is likely to cause this separation. Note that the overall high electrical conductivity masks these lithological differences in the background ERT images. This fine layer is a hydraulic barrier that forces freshwater to flow even more through the preferential flow path provided by the borehole itself and its surrounding gravel pack. Above the fine layer, the plume expands again due to the larger hydraulic conductivity of the coarser sediments.
During the experiment, the water table as well as the electrical conductivity and the temperature of the borehole fluid were measured manually in all five boreholes. The water table rose about 1.5 m in the injection borehole and about 0.2 m in the surrounding four boreholes. The electrical conductivity log of the central borehole before, during and after injection is shown in Fig. 3. It can be observed that during injection (i.e., about 1 h after start of injection), the saltwater in the borehole was pushed up by freshwater. Shortly after injection stopped (5 h after start of injection) the freshwater filled the entire borehole length, whereas it is visible that the saltwater already entered the borehole in the bottom part (at about 16 m depth) and made its way upwards. Therefore, 1 day after the injection experiment, the fluid electrical conductivities in the central borehole were practically back to their initial values, with small differences between 8 and 14 m depth still visible. The electrical conductivities of the fluid in the four corner boreholes showed only small changes that nonetheless indicate that part of the freshwater bulb also reached the outer boreholes.

Synthetic experiment
In order to investigate the behavior of the injected freshwater bulb, and assess in particular the influence of the subsurface hydraulic properties on the bulb evolution, we performed a synthetic study based on the field experiment. This was undertaken using a density-dependent flow and transport simulator. Given the computational burden of the simulations and our goal of examining in detail some of the governing parameters, we did not use a data assimilation approach at this stage, opting instead for analyses of specific scenarios. We considered four scenarios of hydraulic conductivity distribution, and compared the simulated results to each other and with the field evidence in order to gain some first insights on the dynamic response of the hyper-saline-freshwater system.

Flow and transport modeling
For the coupled flow and transport modeling of the freshwater injection experiment, we used a 3-D density-dependent mixed-finite element-finite volume simulator (Mazzia and Putti, 2005). This algorithm was shown to be very effective in the presence of advection-dominated processes or instabilities in the flow field induced by density variations (Mazzia and Putti, 2006). Here, groundwater flow is described by where v is the Darcy flux or velocity, K s is the saturated hydraulic conductivity tensor, ψ is the pressure head and z the elevation head. The hydraulic conductivity is expressed in terms of the intrinsic permeability k and the properties of the fluid as where ρ 0 is the density of freshwater, g the gravitational acceleration and µ 0 the viscosity of freshwater. For densitydependent flow, the density and viscosity of the solution are strongly dependent on the concentration of the solution: Here c is the normalized concentration (i.e., the ratio between the concentration of the solution and the maximum concentration) and and are the density and viscosity ratios, respectively, defined as where ρ s and µ s are the saltwater maximum density and viscosity, respectively. In our case, the density and viscosity ratios are = 0.084 and = 0.28, respectively (see also  Table 1). For the exponential laws in Eq. (3a) and (b), we used a linear approximation (i.e., ρ = ρ 0 (1 + c) and µ = µ 0 (1 + c)) to reduce the computational cost while introducing only a negligible inaccuracy. The mass conservation equations for the coupled flow and transport model can be written as (Gambolati et al., 1999) where S s is the specific storage, t is time, η z is the unit vector in z direction, φ the porosity, q * is a source (positive)/sink (negative) term, v is the Darcy velocity, D is hydrodynamic dispersion, c * is the normalized concentration of salt in the injected/extracted fluid and f is the volumetric rate of injected (positive)/extracted (negative) solute that does not affect the velocity field (Mazzia and Putti, 2006). For the flow and transport model we used a 3-D mesh (Fig. 8) with about 57 000 tetrahedral elements and 10 000 nodes. The size of the mesh was a good compromise between mesh resolution and computational effort. The computational domain extends for 20 m in the x and y directions and 15 m in z direction, starting at 5 m b.g.s., thus representing only the saturated zone. This choice focuses our attention on the processes of interest and reduces dramatically the numerical complexity of modeling coupled flow and transport processes in variably saturated porous media. However, because a water table rise was observed in the boreholes during the injection experiment, we needed to account for this pressure transient in the flow and transport model. Thus, we simulated a comparable injection experiment using a 3-D variably saturated flow simulator (Paniconi and Wood, 1993). The changing pressure values due to the water table rise at 5 m depth were then taken as top boundary conditions for the fully saturated flow and transport model.
In addition to the boundary condition described above for pressure and with c = 0, we set Dirichlet conditions also on the lateral boundaries with a hydrostatic pressure, according to the concentration dependency ψ = − (1 + c) z and Neumann no-flow conditions at the bottom of the mesh. The flow and transport parameter values are given in Table 1. The injection borehole was modeled as a preferential flow path by giving the corresponding cells a large value of hydraulic conductivity. Also the borehole backfill material was included in the simulation by giving it a slightly higher hydraulic conductivity than the surrounding aquifer material. The salt

Injection parameters
Injected volume V mod 20 m 3 Injection duration 3.5 h concentration was given as normalized concentration with a value of 1.0 for the saltwater and 0.0 for the injected freshwater. The initial conditions for the concentration in the aquifer were set to honor the transition zone observed in the borehole fluid conductivity log (Fig. 2). The conditions for the injection were set by giving the cells that represent the injection chamber (between 13 m and 14 m b.g.s.), a pressure head ψ 2 m higher (from 15 to 16 m). To simulate the emptying of the tank, the pressure head decreases over time, calibrated after the measured injection rate in the field.
The immediate increase of the injection rate, observed in the field experiment, was modeled by a "de-clogging" effect of the material closely surrounding the injection chamber (i.e., representing the backfill material). This was done by increasing the hydraulic conductivity of the corresponding cells by about 1 order of magnitude after a correspond- ing time (i.e., about 5000 s). The simulated and true injection rates are compared in Fig. 9.
Dispersive processes play a minor role for the relatively short timescale of the experiment. In fact, several dispersivity values were tested and compared (modeling results not shown here); their influence is not significant over the short timescale considered here. Thus, only advective transport is studied.
To investigate the influence of heterogeneous hydraulic conductivity distributions in the aquifer, four different scenarios were simulated, including one homogeneous model and three different layered models, with a fine (clay-silt) layer between 10.5 and 11.5 m depth ( Table 2). The hydraulic conductivity values for the different scenarios were calibrated manually.

Simulation of ERT monitoring
In order to compare, at least in a semi-quantitative manner, the observed ERT inversions with the results of the synthetic study, it is necessary to convert first the simulated normalized salt concentration from the flow-transport model into bulk electrical conductivity, for example through Archie's law relationship (Archie, 1942), here expressed for saturated sediments: where σ b is the bulk electrical conductivity, a is a tortuosity factor, σ w is the electrical conductivity of the fluid and m is the cementation exponent. The formation factor F = a/φ m accounts for the pore space geometry. Due to the high salinity of the groundwater in the present case, surface conductivity is assumed to be negligible, and thus Archie's law is safely applicable. Since core data were available from one of the boreholes, it was possible to calibrate Archie's law in the laboratory with F = 4.6. The next step is to simulate the field data that would be acquired given the simulated bulk electrical conductivity. For the 3-D electrical forward modeling, we used the same approach as Manoli et al. (2015) and Rossi et al. (2015). The electric potential field, , for a current injection between electrodes at r S+ (current source) and r S− (current sink) is calculated by solving the Poisson equation together with appropriate boundary conditions, where σ b is the given electrical conductivity distribution, I is the injected current strength and δ is the Dirac delta function. The mesh for the geoelectrical modeling includes the unsaturated zone, and the top boundary of the mesh (at z = 0 m) was set as a Neumann no-current boundary condition. For the lateral and bottom boundaries we used Dirichlet boundary conditions. Therefore, the mesh size was expanded in all directions with respect to the hydraulic mesh, so that the influence of the fixed voltage boundary conditions on the current lines was negligible. The final step was to process and invert the synthetic ERT data in the same way as the field data. Upper layer 5 × 10 −5 m s −1 5 × 10 −5 m s −1 1 × 10 −3 m s −1 1 × 10 −3 m s −1 Middle layer 5 × 10 −5 m s −1 1 × 10 −6 m s −1 1 × 10 −6 m s −1 Bottom layer 5 × 10 −5 m s −1 5 × 10 −5 m s −1 1 × 10 −5 m s −1 1 × 10 −5 m s −1

Moment analysis
In order to provide a more quantitative comparison between the field and synthetic experiments, we analyzed 2-D moments as defined for example by Singha and Gorelick (2005): where M ij is the spatial moment of order i, j between 0 and 2, x and z are the Cartesian coordinates and dx and dz the pixel sizes. is the integration domain of interest. The zeroth moment represents the total mass in the system while the vertical first moment, normalized with respect to mass, defines the center of mass in the z direction. The second moments relate to the spread around the center of mass.

Results and discussion
As a first step, let us consider the results of the synthetic study. Figure 10 shows the salt concentration of the flow and transport simulations for scenario 4, which represents the most complex parameterization of the aquifer and is assumed to be most realistic for the test site (see the site stratigraphy reported in Fig. 2). A general upward motion of the injected bulb is visible, with the highest velocities occurring within the injection hole. After some time, the freshwater starts to enter the aquifer along the entire borehole length. Although its density is much less than the density of the surrounding saltwater, the freshwater also moves downwards within the borehole, pushed by the pressure gradients. The 1.2 m thick fine-material layer also plays a clear role in the bulb dynamics. This is expected. In correspondence to this layer, the flow only takes place along the borehole and the backfill material. Above the fine layer the plume expands laterally into the aquifer. Also the transition between the saltwater and the upper freshwater layer above 7.4 m depth moves entirely upwards since the overall movement in the model domain is upwards. One can also observe in the simulation results the tilting of the freshwater-saltwater interface in the lower part of the borehole as well as below the groundwater level, as described by Ward et al. (2007Ward et al. ( , 2008. The higher the ratio of hydraulic conductivity between the two layers, the stronger is the tilting, as predicted by Ward et al. (2008) (results not shown here). Figure 11 shows the inverted images for four different subsurface scenarios at time 4.2 h after start of injection for the flow and transport simulations and the synthetic ERT monitoring (see Table 2 for definition of the scenarios). The figure clearly shows the dramatic influence of the hydraulic conductivity distribution on the shape of the freshwater bulb, both in the "real" images and in the corresponding inverted ERT images. Scenario 4, which includes the fine layer, is closest to the field results as already discussed above. However, scenario 3, with just two layers, shows a similar behavior in terms of plume development. In general, given the strong influence that hydraulic conductivity has on the results, it is conceptually possible to try and infer the site's hydraulic properties on the basis of the freshwater injection experiment. However, it is also apparent that calibrating in detail the true hydraulic conductivity distribution in the field experiment starting from the ERT images alone may be a very challenging task. In fact, while some main features are clearly identifiable, other smaller details may prove difficult to capture.
Indeed, the governing hydraulic effect comes from the different conductivities of the upper and lower parts of the aquifer (scenarios 1 + 2 vs. 3 + 4), and the fine layer does not play such an important role as expected a priori. From the simulation results it is difficult to say whether scenario 3 or scenario 4 is closest to reality. However, for scenarios 1 and 2 ERT clearly overestimates the extent of the freshwater plume, whereas for scenarios 3 and 4 the plume extension is reconstructed quite well, in particular in the deeper region (Fig. 10).
It is instructive to examine in detail (Fig. 12) the similarities and differences between the ERT field data and the reconstructed ERT images from the simulation scenario that visually appears better than the others (scenario 4). The simulated ERT images show the same general behavior in response to the injection process and associated plume development as the ERT field results. In the field ERT images the freshwater body disappears much faster. After 24 h, although in the field ERT images the freshwater bulb is hardly visible, the simulation still shows its presence. It should be noted that in the simulations the boundary condition at the well is imposed as a Dirichlet (head) condition, so flux is computed depending on the applied head. We applied the head as actually measured in the injection tank. Consequently, the flow is never zero, not even at the end of the experiment. On the other hand, the tilting of the freshwater-saltwater interface Figure 10. Flow and transport modeling results at different times (in hours after start of injection) for scenario 4 (see Table 2). as seen in the flow and transport model results is much less visible in the ERT images.
The imaged resistivity changes in the field experiment show less contrast than in the synthetic study. The salinity difference between the freshwater and the saltwater is very large and thus so is the NaCl concentration. Within this range, the electrical conductivity of the water might no longer follow a linear relation with concentration (e.g., Wag-ner et al., 2013), while here it is assumed to be linear. This can lead to a shifting in the contrast when the concentration is converted into electrical conductivity.
Note also that the gradual change of electrical conductivity in the transition zone (i.e., between 5 and 7.4 m depth) is not visible in the ERT images (Fig. 11). In the transport simulations it can be seen that this zone also moves upwards in the aquifer and becomes thinner (Fig. 10). Figure 12. Results of synthetic ERT experiment for selected times (in hours after start of injection) for scenario 4 (see Table 2). Black diamonds denote the position of electrodes.
Another difference between the field and the synthetic ERT results is the sharpness of the freshwater body; the boundaries appear smoother in the field study. Although dispersion effects were not further investigated in this study, a higher value of α L and α V in the simulations would obviously lead to a smoother gradient across the plume boundaries. On the other hand, in the field results this may also be partly explained by the fact that one ERT measurement frame took about 40 min, and since the overall plume migration was relatively fast, the process is to some degree smeared in the inverted images. Figure 13 shows the spatial moments (0th moment: total mass; 1st moment: center of mass) of the freshwater bulb for the field and synthetic ERT inversion results, as well as the "true" moments from the flow and transport model (see Sect. 3.3). The total mass is well recovered by the synthetic ERT results (using backwards the same Archie's law parameterization used in the forward modeling). However, the field ERT underestimates the total mass. While this is a known characteristic of moment analysis applied to ERT data for tracer tests (e.g., Singha and Gorelick, 2005), in this specific case it looks likely that the chosen Archie's law parameters are not fully adequate to represent the electrical conductivity-salinity relationship. Considering that even linearity of Ohm's law is questionable at the high salt concentrations observed at the site, one could also question the overall validity of Archie's law. Note that all other factors normally contributing to bad ERT mass recovery under field conditions are the same in the synthetic and the true case, and thus cannot be called into play. Figure 13. Spatial moments for the field ERT data, synthetic ERT data, and the true data from the flow and transport model. The moments for the true field were calculated in 3-D while those for the ERT tomograms were calculated in 2-D. The field ERT data are separated into the two borehole planes. (a) shows the total mass in the system, normalized, and (b) is the center of mass in the vertical direction.
In contrast to the total mass, the vertical center of mass is, despite some early oscillations, well recovered also for the field data. This, however, is known to be a very robust indicator (e.g., Binley et al., 2002;Deiana et al., 2007Deiana et al., , 2008.
Overall, and in spite of the differences described above, the comparison between observed and modeled ERT images is satisfactory, particularly in the face of uncertainties concerning the heterogeneities of the real system that could not be investigated in extreme detail. In addition, we cannot exclude the possibility that the linearity of the current flow equation may be violated in such a highly conductive envi-ronment, thus leading to inconsistencies between field reality and theoretical assumptions.
Despite the above limitations, the comparison shows that ERT imaging is a viable tool for monitoring freshwater injection in a hyper-saline aquifer. This, by itself, was not an obvious result. The ERT dataset was collected under extreme, challenging conditions. Even so, the ERT data are of fairly good quality considering that we retained only data that passed a fairly strict reciprocity check, knowing that larger reciprocity errors are likely to be related to nonlinear current effects occurring in such high electrical conductivity environments. The study also indicates how an accurate coupled model can mimic in an effective manner the behavior of the observed freshwater bulb that was injected into the domain, and this too was not self-evident.

Conclusions
In this paper we present a hydro-geophysical approach that can be used to study freshwater injections in saline aquifers. In particular the approach is used to monitor and describe a freshwater injection experiment conducted in a hyper-saline aquifer in the Molentargius Saline Regional Park in the south of Sardinia (Italy). The experiment was monitored using time-lapse ERT in five boreholes. A numerical study of the experiment (density-dependent flow and transport modeling in conjunction with ERT simulations) was carried out to investigate the plume migration dynamics and the influence of different hydraulic conductivity parameterizations. The numerical algorithm of the coupled flow and transport model proved to be stable and accurate despite the challenging conditions.
The results demonstrate the feasibility and benefit of using a combination of (a) time-lapse cross-borehole ERT and (b) numerical modeling of coupled flow and transport to predict the same ERT results. The comparison between measured and simulated ERT images was used as the key diagnostics aimed at estimating the system's governing parameters and consequently describing the saltwater-freshwater dynamics. More sophisticated data assimilation techniques can be used to further refine the presented approach in future work. We can conclude from the present study the following: a. The complex dynamics of hyper-saline-freshwater systems can be tracked using high-resolution spatially extensive time-lapse non-invasive monitoring. On the contrary, traditional monitoring techniques alone (e.g., conductivity logs, as in Fig. 3) give only a very partial image, largely inconclusive to understand the system dynamics.
b. Numerical modeling of these coupled systems is very challenging due to the presence of strong density/viscosity contrasts and large hydraulic conductivity heterogeneities. The latter, in particular, largely control the dynamics of the saltwater-freshwater interaction. In absence of a robust numerical model, it is impossible to estimate the impact of hydraulic heterogeneity on this dynamics. c. A detailed comparison between field data (here, ERT time-lapse images) and modeled data of the same type enables a better understanding of the behavior of a freshwater bulb injected into a hyper-saline environment.
Our study also serves to highlight some of the weaknesses that should be addressed in future work: -Fine-tuning of geophysical constitutive relationships, hydraulic and transport parameters, and system heterogeneities needs to be improved. We managed to bring the match between field and synthetic data to an acceptable level with relatively small effort, but it is very difficult to improve the match further. For instance, in the case presented here the injected freshwater bulb "disappears" from the real ERT images faster than in the simulation results. Also, the mass balance is honored easily in the simulations, whereas in the real data lack of mass is apparent. All of this points towards a number of aspects that could be improved in the data matching. However, the target parameters to be modified for this improvement are not easy to identify, given their very high number and complex nature. Among these, there are hydraulic parameters and dispersivities, and their spatial heterogeneities, as well as also Archie's law parameters. This task is likely to be challenging even in a rigorous data assimilation framework, and equifinality of model parameterizations is likely.
-The extreme hyper-saline system considered here is likely to exceed the limits of linear relationships between current and voltage (Ohm's law) as well as between electrical conductivity and salinity. Therefore, a full nonlinear analysis should be conducted, particularly concerning the electrical behavior of the system. In absence of this, we have to limit ourselves to a semiquantitative interpretation, as shown here.
Finally, with regards to practical aspects of freshwater injection and monitoring in saline aquifers, we can draw the following conclusions: -Although in typical ASR applications the contrasts of density and salinity are usually smaller, this study shows that time-lapse ERT is a powerful monitoring tool for this (and also other) type of hyper-saline applications.
ERT can provide spatial information that is unattainable using traditional monitoring techniques (e.g., in boreholes).
-The movement and mixing of the freshwater plume can be very fast; thus, any ERT monitoring must adopt con-figurations for quick measurements (e.g., in the conditions represented in this study an acquisition time of less than 30 min is recommended).
-In hyper-saline systems, measuring reciprocity may not be the ideal error indicator since nonlinear phenomena may be triggered, or during the time between the normal and reciprocal measurement the system may have already changed, thus invalidating the reciprocity check.
The example shown in this paper shows how the joint use of ERT imaging and gravity-dependent flow and transport modeling give fundamental information for this type of study.

Data availability
Measured raw cross-borehole time-lapse ERT field data, additional field data, inverted ERT field data as well as the modeling data in terms of the concentration distribution of the density-dependent flow and transport model and the inverted synthetic ERT monitoring results can be accessed at doi:10.5281/zenodo.322630 (Haaken et al., 2017).