Experimental observations of aquifer storage and recovery in brackish aquifers using multiple partially penetrating wells

Aquifer storage and recovery systems using multiple partially penetrating wells (MPPW-ASR) can form a viable solution to the problem of freshwater buoyancy when using brackish aquifers for freshwater storage. This study presents the result of a series of laboratory experiments that aimed at visualizing the shape of freshwater bodies injected into a brackish aquifer and determining the effect on the recovery efficiency (RE) of several MPPW-ASR operational variables. A model aquifer was built in a Plexiglas tank using glass beads and water was injected and abstracted through point and vertical wells, which were operated in various combinations. Numerical models were used to support the interpretation of the time-lapse photographs, and showed that three-dimensional flow effects had to be considered for a correct interpretation of the visible dye patterns. Upward migration of both fresh (during injection) and brackish water (during recovery) along the vertical wells was observed, indicating that the role of well infrastructure as conduits is a critical design criterion for real-world systems. Gravitational instabilities formed when freshwater did not extend all the way to the top of the aquifer, and this negatively impacted the RE by causing greater mixing. The positive freshwater buoyancy led to freshwater bodies that became narrower with depth, and the formation of thin, elongated buffer zones along the aquifer top in multicycle experiments. Up-coning below abstraction wells resulted in lower RE values, reinforcing the potential of scavenger wells to enhance MPPW-ASR system performance.


Introduction
Managed aquifer recharge refers to the broad suite of methods to maintain, enhance and secure groundwater systems, as well as to protect and improve water quality (Dillon et al. 2019). Aquifer storage and recovery (ASR) is one form of managed recharge that involves the injection and recovery of water into an aquifer using a single well. Many successful ASR operations around the world are testimony to the feasibility of this technology for the underground storage of water during times of surplus, followed by the recovery to meet the demand shortfall during dry periods or emergency events. Over 175 ASR wellfields are operational in the USA (Dillon et al. 2019), followed by Australia with 35 operational sites in Adelaide alone in 2016 (Kretschmer 2017). In Europe 12 well injection sites were reported operational in 2013 (Sprenger et al. 2017) and new sites have been implemented since . The technology is also increasingly implemented elsewhere such as on the Arabian Peninsula ) and in Latin America (Bonilla Valverde et al. 2018).
Even though ASR practices in saline aquifers can be traced back to as early as around AD 600 in India (Sakthivadivel 2007), storing freshwater in saline groundwater is associated with a number of problems that have thus far hampered the widespread application of ASR in coastal areas. The two main difficulties are (1) the mixing of the injected freshwater with ambient saltwater and (2) the floating to the top of the aquifer of the freshwater due to its lower density compared to the saltwater (Bakker 2010;Esmail and Kimbler 1967;van Ginkel et al. 2014;Zuurbier et al. 2013).
Various studies have investigated the influence of these processes on the recovery efficiency, i.e. the volume of usable water that can be abstracted as a fraction of the volume that was injected. Esmail and Kimbler (1967) conducted a series of laboratory experiments and found that density effects impacted more negatively on the recovery efficiency than did mixing effects. Furthermore, wider mixing zones were found to suppress the rate of tilting of the boundary between the dense and less-dense fluids used in their study. While Esmail and Kimbler (1967) used separate laboratory setups to study the effects of mixing and density individually, a follow-up study by Kumar and Kimbler (1970) used a single, pie-shaped model representing a 45°sector of a circle. The experimental results were used to validate an analytical solute transport model of the flow of miscible fluids for a radial system, which was then used to determine the effect of various parameters on the recovery efficiency. It was found that the recovery efficiency increased with the number of injection-abstraction cycles. Conversely, the recovery efficiency decreased with increasing density difference, aquifer permeability or thickness, and duration of the storage phase.
In a series of papers, Ward et al. (2009Ward et al. ( ), (2008, (2007) explored the influence of the controlling variables on ASR in brackish aquifers using numerical models. Their work showed that a low pumping rate, high hydraulic conductivity, low dispersion, and long storage phase duration all increased the importance of density effects. While dispersive mixing in principle had a negative impact on the recovery efficiency, Ward et al. (2007) concluded that a widening of the mixing zone reduced density gradient and thus the negative impact of the density effect. This confirms earlier assertions by Esmail and Kimbler (1967) and explains Kumar and Kimbler's (1970) finding that the recovery efficiency increased with increasing dispersion coefficient. Ward et al. (2008) further found that a greater difference between high and low permeability in a layered aquifer and a high vertical anisotropy in a homogeneous system suppress the vertical flow due to the density effect and thereby increases the recovery efficiency. Bakker (2010) used a simplified modelling framework (i.e. based on a sharp interface approach that neglects dispersive mixing), which allows for a first-order assessment of the recovery efficiency based on the ratio of the well discharge and the product of the hydraulic conductivity, the square of the aquifer thickness, and the dimensionless density difference.
Several data have been collected during field trials in, for example, Virginia (Cederstrom 1947) and Texas (Moulder and Frazor 1957) in the USA, and more recently in the Netherlands (Zuurbier et al. 2014). The latter system used multiple partially penetrating wells (MPPW-ASR), which enabled the injection of freshwater into the deeper part of the aquifer, while abstraction remained restricted to the wells screened in the upper part of the aquifer. Comparative numerical model simulations showed that this system can achieve recovery efficiencies of over 40% vs~15% for a conventional fully penetrating well. Modelling showed that almost 60% of the injected water in a year could be recovered by the MPPW in additional cycles (Zuurbier et al. 2014).
To date, MPPW-ASR systems in saline aquifers have not been studied as extensively as single well ASR systems. A need remains to investigate how operational design factors like well arrangement, pumping scheme and injection rate can be optimised to increase the recovery efficiency. In the present study, a series of sand tank laboratory experiments was conducted to visualise the injected freshwater during the injection, storage and recovery phases. The aim was to investigate how the aforementioned operational design factors affect the recovery efficiency. Numerical models were used to assist in the interpretation of the visual images of the dye patterns.

Tank setup
For the laboratory experiments, a model aquifer (hereafter referred to simply as aquifer) was created in a 1.87-m-long, 0.05-m-wide and 0.546-m-high Plexiglas tank (Fig. 1). Two stainless steel grids covered with mesh fabric were built into the tank to section off 0.035 m on each side, which allowed having two water reservoirs for controlling the heads in the aquifer. The aquifer was made up by spherical glass beads with a diameter of 1-1.3 mm. The tank was filled up to a height of 0.275 m by pouring the glass beads into the tank and levelling them to create a horizontal surface. The water reservoirs on both sides of the tank were then slowly filled with distilled water to saturate the glass beads from below, in order to avoid the formation of air bubbles. To achieve confined conditions, the aquifer was covered by two Plexiglas panes (which had been cut to fit snug into the tank) that were inserted after the beads were fully saturated and the water level was at the same height as the top surface of the aquifer. The panes were sealed using a commercially available sealing compound (plastic-fermit).
The horizontal hydraulic conductivity (K) of the glass beads was determined from a series of horizontal flowthrough experiments. A total of nine different head gradients (I) across the tank was considered by varying the water level in the left reservoir. Once it was verified that the flow rate (Q) had become steady, it was determined by measuring the volume of water flowing from an outlet on the right during a time span of 2 min. Using Darcy's Law, the hydraulic conductivity was calculated as K = Q/(A × I), in which A is the crosssectional area of the outflow face (A = 0.014 m 2 ). From the nine experiments, using 0.004 < I < 0.019, the hydraulic conductivity was found to be K = 0.93 m/min with a standard deviation of 0.04 m/min.
During the process of levelling the top surface of the glass beads, as well as during the installation of the confined layer, some compaction took place. Hence, ex-situ measurements of porosity are inaccurate and the porosity was determined by treating it as a calibration parameter in the numerical model of one of the laboratory experiments (see below for further details). Ex-situ porosity measurements of the glass beads yielded a range of values from 0.34 to 0.41, depending on the degree of compaction achieved by shaking the measuring cylinder. The theoretical range for monosized spheres lies in the range 0.36-0.4 (Dias et al. 2004).
Three different MPPW-ASR well configurations were assessed: (1) two vertical wells, (2) four point wells, and (3) two point wells. The two vertical wells consisted of Plexiglas tubes that were inserted into the aquifer from the top through two holes in the confining pane. Well screens were created by inserting eight perforations into the tubes over a length of 0.105 m (i.e., a 0.015-m spacing between the perforations). The wells were inserted to different depths such that the bottom of the shallowest screen was 0.015 m above the top of the deepest screen (Fig. 1b). The point-wells were created using Plexiglas tubes, inserted 0.025 m into the aquifer from the back of the tank. The point-wells were filled with glass wool to prevent glass beads from entering the tubes.
To abstract and inject water, a peristaltic pump (ISMATEC, Reglo ICC) with four individually controllable channels was used. Due to the confined conditions in the aquifer, the pumps were able to extract water more easily than to inject it. Because this made it difficult to inject at an exact, predefined rate, the actual injection (and extraction) rates were determined by measuring the injected and extracted water volumes over time. The water levels in the reservoirs were kept at a constant height of 0.30 m using Mariotte's bottles when water was being abstracted from the well system. During injection phases, the reservoirs' water level was maintained using a pump that skimmed water from the top of the water columns through an inlet tube that was placed at 0.30 m above the base of the tank (Fig. 1a). By setting the pump to a pump rate much higher than the freshwater injection rate, it was ensured that the reservoir water level never rose above the height of the inlet tube opening.
Depending on the experiment, the glass beads were saturated with purified (electrical conductivity <10 μS/cm) or brackish water. The latter was created by dissolving 6.00 ± 0.3 g/L sodium chloride in purified water that had been heated to 80°C. The brackish water was tinted orange using sulforhodamine G. The freshwater that was injected through the lower screens was dyed yellow using uranine AP, while the freshwater injected through the upper screens was coloured blue using indigo carmine. All dyes were added at a concentration of 0.08 g/L. The experiments were captured by a wide-angle camera (GoPro HERO5), which took a picture every 10 s. The camera was mounted on a tripod that was placed on a fixed position in front of the tank. Image postprocessing was done to correct for minor camera tilt and In operational ASR management, the recovery efficiency (RE) is defined based on one or multiple water quality threshold criteria of the extracted water (Kimbler et al. 1973;Merritt 1986). For the laboratory tank experiments, a threshold value for the specific conductance (SC, i.e. the electrical conductivity at 25°C) of 2.75 mS/cm was used. This value was chosen based on the threshold value of SC = 2.790 mS/cm as specified in the German drinking water ordinance (BGBl 2016). During recovery cycles, the electrical conductivity of the abstracted water was measured at the outlet of the sand tank using flow-through PendoTECH Single Use Conductivity Sensors (accuracy of 5%). An automated temperature correction was performed to convert the electrical conductivity to SC. Only two channels could be recorded simultaneously, which necessitated exchanging the sensors connected to the conductivity meter during the course of an experiment. The brackish water had an SC = 10.75 ± 0.5 mS/cm, corresponding to a salinity of S = 6.08 ± 0.3 g/L (Clesceri et al. 1999). The distilled water had a density of 997.5 kg/m 3 , while the density of the brackish water was 1,001.8 kg/m 3 . The density was measured using a portable density meter (Anton Paar, DMA35). The experiments were all conducted at room temperature, which ranged between 21 and 23°C.

Experiments
The three main operational variables of a MPPW-ASR scheme that were varied between experiments were the (1) well arrangement, (2) pumping rate, and (3) pumping scheme. Various combinations of these variables led to 16 experiments, which are summarized in Table 1. Experiments 1, 2 and 3 were carried out using distilled water as the ambient and injected fluid. This was done to exclude the confounding effects of density-driven flow on the flow field. All other experiments were conducted using brackish water as the ambient fluid.
The three different well arrangements (vertical, two-point and four-point) were described in the preceding. Two injection/recovery rates were used, which differed by a factor of two: In all experiments approximately 3 L of water was injected (Table 1) either over a time span of 60 or 120 min. For the experiments with the 120-min injection cycles, the recovery rates were then roughly half (with differences being attributable to the accuracy to which the pump rates could be controlled) of their 60-min equivalents. In what follows, experiments with the higher injection/abstraction rates will be indicated as category H, those with the lower rates as category L experiments, respectively.
Two different pumping schemes were applied. For pumping scheme 1, infiltration was only through the deep screen if the vertical (V2, see Fig. 1b) or two point wells arrangement (screen P4) was used, or through the two deepest screens in the case of the four-point well arrangement (screens P3 and P4). Recovery was only through the shallow screen of the vertical and two point wells (V1 and P1, respectively), and the two topmost screens of the four point wells (P1 and P2). In pumping scheme 2, both injection and abstraction were through the deep and the top screens of the vertical (V1 and V2) and two point wells (P1 and P4), or through all the screen of the four point wells (P1-P4). The injection rates increased with depth to create as wide a freshwater body as possible. Conversely, during recovery, the abstraction rates were highest near the top of the aquifer. The distribution of the injection and extraction rates across the individual well screens for the pumping schemes is provided in Table S1 of the electronic supplementary materials (ESM).
The majority of the experiments were carried out with all the wells being active over the entire experiment duration, regardless of whether they were abstracting fresh or brackish water. As will be shown below, the brackish water reached the deepest wells first and once they were abstracting brackish water, they were acting as scavenger or interception wells. Scavenger wells are wells that sit below the injected freshwater and remove upward flowing saltwater with the aim to slow down the salinization of the freshwater abstraction wells (Saravanan et al. 2014;van Ginkel et al. 2014;). To simulate a MPPW-ASR scheme without scavenger wells, experiment 14 (which was selected because it had the highest RE) was repeated (experiment 15), but rather than all wells pumping continuously, individual wells were deactivated once the threshold value of SC = 2.75 mS/cm (corresponding to a salinity of 1.44 g/L) was exceeded.
In addition to the single-cycle experiments in Table 1, three experiments were conducted to evaluate the effect of multiple infiltration, storage and extraction cycles on the RE. Experiments 14 and 6 were selected as they had the highest and lowest RE, respectively, while experiment 15 (as 14 but without scavenger wells) was also conducted as a multicycle experiment. The multicycle variants are denoted by the experiment number with the letter M as a suffix. Five consecutive injection-storage-recovery cycles were considered. Since the experiments consequently lasted several hours, the storage cycle duration was decreased to 5 min in order to keep their total duration within practical limits. The details of the experiments can be found in Table S2 in the ESM.

Numerical modelling
Numerical simulations were conducted using SEAWAT (Langevin et al. 2008) using FloPy (Bakker et al. 2016) as the model interface. The three-dimensional (3D) model domain had the same dimensions as the volume occupied by the glass beads in the tank. The number of layers, rows and columns was 55, 5 and 90, respectively, corresponding to a cell dimension (height × width × length) of 0.005 × 0.01 × 0.02 m. No-flow and zero-concentration gradient conditions were specified for the top, bottom, front and back boundaries. Along the left and right boundaries, the heads were specified at 0.3 m, in accordance with the setup of the laboratory experiments. The same head value was used as the initial head across the entire model domain. The concentration boundary condition on the left and right boundaries depended on the flow direction, with inflowing water having a concentration of 6.08 g/L and outflowing water having a concentration equal to that of the boundary cell. The initial salt concentration was also 6.08 mg/L.
The two vertical wells were represented by 16 model cells, the locations of which coincided with the points where the vertical tubes had been perforated, or where the point wells were inserted through the back wall of the tank. The duration of the model simulation and the injection, storage and recovery cycles was varied depending on the experiment designs listed in Tables S1 and S2 of the ESM. One-minute stress periods subdivided into six flow time steps were used.
The experimentally determined hydraulic conductivity of K = 0.93 m/min was used. It was assumed that the glass beads were isotropic, except where the vertical wells were located. The cells intersecting these wells were assigned a vertical hydraulic conductivity of K = 9.3 m/min. This was done to account for preferential vertical spreading of the dye tracer that was observable in some experiments, as will be shown below. This hydraulic conductivity value was deemed appropriate after running a number of trial simulations. Similarly, the hydraulic conductivity (horizontal and vertical) of the top model layer was increased by a factor of two compared to the experiment value to account for the voids between the top surface of the glass beads and the confining pane, which also caused observable preferential spreading of the injected water. The remaining model parameters are listed in Table 2.
Experiments 3 and 14 were chosen for the calibration of the numerical model. First, the porosity was adjusted by comparing the dye pattern to the modelled extension of the injected water for experiment 3. The resemblance between the model and the observed dye patterns was assessed by overlaying the  simulated isoconcentration lines over the photographs of the tank. Second, the dispersivity was adjusted to optimize the match between the modelled and the measured salinity versus time measurements for each of the abstraction wells in experiment 14. In order to compare the laboratory-measured specific conductances to the numerical results, the measurements were converted to salinities using the equations in Clesceri et al. (1999). The choice of solver for the advection term of the transport equation had a significant effect on the model outcomes. The results using the finite-difference solver displayed too much numerical dispersion and the method-of-characteristics solvers showed noisy breakthrough curves that precluded accurate comparison with the measured concentration-time curves. The only solver that produced satisfactory results was the TVD solver, which was used for all simulations. The TVD solver has been shown to perform adequately when simulating buoyant freshwater in laboratory experiments (Goswami et al. 2012).
Particle tracking calculations were used as an additional aid to visualize the flow of water in the aquifer and were done using MODPATH Version 6 (Pollock 2012). Particles were released with 1-min intervals from four locations within each cell containing a well during the injection stage. The particles were released at a distance of 0.02 from the front face of the model. Visualization of the 3D flow field in combination with the photographs showing the spreading of the dye tracers is inherently difficult. Therefore, only the particles that were within 0.01 cm (i.e., one model cell width) are shown in the figures that follow, meaning that not all particles are visible in the figures. This approach was chosen to provide an indication of the advective transport pathways of the injected dye tracer. It should be borne in mind that the laboratory visualization experiments are also influenced by diffusive and dispersive spreading of the dye tracer.

Flow dynamics and salinity patterns
The comparison between the numerically simulated and laboratory-measured results for the calibration experiments (experiments 3 and 14) is shown in Figs. 2, 3 and 5. The full time-lapse sequences of the experiments are included as video files (ESM2-ESM8). These also show detailed information about the injection and abstraction points, rates and duration.
The lateral spreading of the dyed tracer in experiment 3 ( Fig. 2a-d) was matched using a porosity of n = 0.38, which lies well in the range of 0.34-0.41 that was found from volumetric measurements in cylinders. The match between the simulated salinity of the abstracted water with the measurement data (Fig. 3) was best for a longitudinal dispersivity of α L = 7 × 10 −4 m. A reduction of this value did not cause steeper concentration versus time curves. The value is at the lower end of the range of 5 × 10 −4 < α L < 3 × 10 −3 m reported by Jakovovic et al. (2011) for homogeneous sand (median grain size diameter of 4.6 × 10 −4 m). This result seems reasonable, as the higher sphericity of the glass beads compared to the quartz sands used by Jakovovic et al. (2011) may be expected to result in less tortuous flow paths, thus suppressing hydrodynamic dispersion.
The shape of the freshwater body was not perfectly symmetrical around the injection well (Fig. 2a-d). Lateral flow effects, due to inherent difficulties of maintaining equal water levels on the left and right of the tank (Werner et al. 2009), are the most likely cause for this. Trial simulations showed that these can already become apparent in the visualised dye patterns for head differences between the two reservoirs on either side of the tank as small as 0.5 mm, which are practically impossible to prevent.
For experiment 14, there was an overall good agreement between the dye tracer patterns and the simulated isocontour lines (Fig. 3e-h). The boundary between the yellow and the orange dye displayed a pronounced inward curvature in the upper part of the tank, just below the region occupied by the blue dye. Trial runs during model calibration showed that the model could replicate this shape at the calibrated porosity value of n = 0.38, but not at lower porosity values. This result indicates that visible dye pattern must have been due to the three-dimensional spreading of the injected freshwater, as illustrated in Fig. 4. After 5 min since the start of the injection, small freshwater bubbles have spread concentrically around the injection points (Fig. 4a, yellow isoconcentration surface). After 10 min (Fig. 4b) the freshwater in the deeper part has already reached the tank wall, but in the upper half of the tank, it has not. It does so after 15 min (Fig. 4c), and after 30 min (Fig. 4d), the pattern becomes largely two-dimensional (2D), although a 3D shape remains recognisable along the aquifer top. This explains the inward curvature of the boundary between the orange and yellow dyes seen at the front of the tank (Fig. 3e-h). This result is not obtained when too low a porosity value is chosen in the model, since the available pore space then becomes so small that the injected freshwater hits the front face much earlier, resulting in almost pure 2D concentration distributions.
The bulging outline of the freshwater in experiment 3 ( Fig.  2a-d) was the result of the injection occurring at four discrete points, and the widening towards the aquifer base was due to the injection rate increasing with depth ( Table 1). The higher injection in the deeper wells induced a vertical flow component near the centre of the aquifer, which is clearly visible from the modelled position of the blue particles. Water that left the uppermost well screens first moved vertically up, as well as horizontally, with the vertical flow component becoming less significant as the water moved further away from the point of injection. The flow away from the two deepest point of injection was mainly horizontal.
The effect of buoyancy on the shape of the freshwater body can be inferred by comparing experiments 3 and 14. Notwithstanding minor differences in the injected volumes (Table 1), both experiments were the same except for the ambient fluid, which was distilled water in experiment 3 and brackish water in experiment 14. The experiments were performed using the 4-point well configuration, using all wells during both the injection and the recovery phase. The vertical flow component in experiment 14 (Fig. 2e-h) was more dominant than in experiment 3. Moreover, as particles move away from the well, they started to move downward, indicating a rotational movement in the transition zone between the freshwater and the brackish water. The horizontal spreading of the injected water was preferentially along the top of the aquifer, whereas spreading in the lower half of the aquifer was much more limited compared to experiment 3. Figure 5 shows the contraction of the injected water volume during the recovery phase. Even though the results are  The measured data are shown as points, with vertical bars representing the error of the salinity measurement. The model outcomes are shown as lines. The different coloured lines represent the four abstraction wells (Fig.  1b). The green area indicates where the salinity is below the threshold value of S = 1.44 g/L (SC = 2.75 mS/cm) shown side-by-side, a direct comparison during the recovery phase is not meaningful as in experiment 3 the extraction rate was highest for the deepest wells, whereas for experiment 14 it was the opposite. In experiment 3, the body of injected water therefore largely keeps it shape, and the time lapse sequence of the recovery phase resembles the injection phase but then in reverse . In experiment 14 on the other hand, the freshwater body changed shape completely as it strongly thinned in the vertical direction. This thinning was caused by the increasing abstraction rate towards the top of the aquifer, aided by the positive buoyancy of the freshwater. As can be seen in the time-lapse video of this experiment, the buoyancy effect also caused a minor but noticeable upward drift of the freshwater during the 30-min storage phase.

Influence of pumping scheme and rate
The effect of the pumping scheme can be evaluated by comparing experiments 6 (pumping scheme 1, injection only through the deepest point-well screen) and 12 (pumping scheme 2, injection through the lowest as well as the shallowest point screen). Figure 6a-d shows the development of the freshwater body during the injection phase of experiment 6. After 15 min an irregular tube-shaped zone of blue dye, emanating just to the right of the tank centre from the newly formed freshwater body, became visible. This represents upward flowing freshwater, which is attributed to the presence of the vertical well at this location. The upward buoyancy of the freshwater and the high permeability represented by the well tube gave rise to this rapid upward movement. At later times, convective salt fingers appeared. In experiment 12 (Fig. 6e-h), the shape of the freshwater body that formed was very similar to that in experiment 6, but the saltwater fingering in the upper part of the injected freshwater was less pronounced.
During the course of the recovery phase of experiment 6, a wide transition zone formed at the base of the triangle-shaped freshwater volume as it started to thin (Fig. 7b) and the effect of the abstraction well was clearly visible from the up-coning of brackish water right below it (Fig. 7c,d). This widening of the transition zone was also observed in the sand tank upconing experiments by Jakovovic et al. (2011) andWerner et al. (2009), which they attributed to enhanced mixing by longitudinal dispersion when the flow and the transition zone movement are in the same direction. Observed patterns were similar in experiment 12 (Fig. 7e-h), but because the lower point well was also abstracting water, the thinning of the lens and the up-coning pattern were less pronounced than in experiment 6. In fact, the faint trace of yellow dye that pointed vertical down from the main freshwater lens, shows that a mixture of fresh and ambient brackish water was being captured by the deep point well.
The effect of the injection and abstraction rates can be inferred from a comparison between experiments 12 and 13 (Figs. 6 and 7), the latter having half the pumping rates as the former. The injected freshwater in experiment 13 (Fig. 6i-l) had a much narrower shape than in experiment 12 (Fig. 6e-h), and also it did not reach all the way down to the aquifer bottom Fig. 4 Three-dimensional representation of modelled isoconcentration surfaces representing 5% (yellow) and 95% (orange) of the ambient water salinity) for experiment 14 after a 5, b 10, c 15 and d 30 min since injection started. The black straight lines represent the model boundaries and the dots indicate the well locations. The viewing direction is for an observer standing in front of the tank towards the left side, viewing towards the centre of the tank's front face in experiment 13 as it did in experiment 12. In both experiments, gravitational instabilities formed in the brackish water that separated the deep and shallow freshwater. During the recovery phase, a horizontally-elongated freshwater body formed in experiment 13 that was much thinner than in experiment 12 (compare Fig. 7i-l and Fig. 7e-h). Figure 8 summarises the experimental outcomes in terms of their recovery efficiencies. The RE was consistently higher for pumping scheme 2 than for pumping scheme 1. Moreover, the RE depended critically on the duration of the recovery phase: The RE was found to be up to 2.1 times higher for the experiments with the higher pumping rates (category H, with H having pumping rates roughly twice as high as L). Also, for category H experiments, the 4-point well system outperformed both the vertical-and the 2-point well system, but for category L experiments, the differences in RE between the well configurations became less apparent. In experiment 15 (not shown in Fig. 8) abstraction was stopped once the threshold SC of 2.75 mS/cm was exceeded, which, as expected, caused a lowering of the RE compared to experiment 14 (56% vs 69%). Figure 9 shows the REs for the laboratory experiments with multiple cycles. As expected, the RE increased with repeated injections. For experiment 14 M, 15 M and 6 M, the greatest increase in RE occurred after the first and second cycle, after which a marginally-increasing or a stable RE was observed in 14 M and 15 M, while in experiment 6 M, the RE increase was still notable after the third cycle. Experiment 15 M, carried out without scavenger wells, proved to have the lowest overall recovery efficiencies. The higher RE after the first cycle in all multiple cycle experiments compared to the single cycle The numerical model approximated the laboratorymeasured RE to a good degree of accuracy (Fig. 10). The reason that most data points plot below the 1:1 line is that the modelled salinity at the abstraction points exceeded the threshold value sooner than during the tank experiment, yielding a lower modelled RE than the observed RE. This is due to the overestimation in the numerical model of dispersive mixing, causing earlier arrival of the salinity front at the wells than in the tank. Even though the TVD scheme showed the least numerical dispersion of all the advective solvers in SEAWAT, a further improvement could not be achieved. A finer grid resolution might abate this issue but the results were deemed satisfactory given the impracticalities associated with longer model run and post-processing times. Whilst the RE of experiments 1, 2 and 3 could not be measured in the laboratory as both the injected and ambient fluids consisted of distilled water, the RE calculated based on the numerical models suggested a RE of approximately 89% for all three experiments.

Recovery efficiencies
For experiments in which buoyancy effects caused strong upward flow of freshwater with salt fingering (like in experiments 6, 12 and 13 described previously, but also experiment 9, which is not shown here), the RE was higher in the models than in the laboratory experiments. The reason lies in the fact that the observed gravitational instabilities were not reproduced by the numerical simulations. Even though the models captured the overall shape and the wide mixing zones that were observed, the loss of freshwater by enhanced mixing due to the fingering flow appears to be underestimated in the numerical simulations.

Discussion
The patterns visualised during the experiments provide direct insight into the processes that influence the storage of freshwater in brackish aquifers as well as the effect of operational factors of MPPW-ASR systems on the RE. The numerical models were able to replicate the main features of all the experiments conducted with a high degree of accuracy. This provided more confidence in the interpretation of the dye patterns and the salinity measurements, but conversely, it validates the use of variable-density flow and transport codes to The results for experiment 12 for the corresponding points in time are also shown (e-h). The results for experiment 13 are shown in the right pane but because the injection rate was twice as small as in experiment 12, the snapshots displayed are after i 30, j 60, k 90 and l 120 min. See the caption of Fig. 2 for a detailed description of the lines and dots predict the spreading of injected plumes in ASR systems. Three-dimensional flow effects inside the tank had to be considered for a truthful numerical simulation of the dye patterns. These could be replicated without the need to consider confounding factors that have been found to be of potential importance in other studies such as boundary wall effects or dye sorption to the glass beads, which can obfuscate the visualization of the flow pattern (Jakovovic et al. 2011), suggesting that these played no significant role in the present study. The results for experiment 12 for the corresponding points in time are also shown (e-h). The results for experiment 13 are shown in the right pane but because the recovery rate was twice as small as in experiment 12, the snapshots displayed are after i 30, j 60, k 90 and l 120 min. See the caption of Fig. 2 for a detailed description of the lines and dots Each of the three groups is subdivided into high-and low abstraction/injection rate (darker vs lighter background colour, respectively). For each subgroup, the darker bar is for pumping scheme 1, the lighter bar is for pumping scheme 2 The buoyancy effect caused by the density difference between the injected freshwater and the ambient brackish water resulted in a strong upward movement of the freshwater, followed by a horizontal spreading of the freshwater along the top of the aquifer (Fig. 2). During the injection and storage phases, a rotational movement of the interface between the two water types was evident and confirmed by the calculated pathlines in the particle tracking simulations. As has been found in other studies (Bakker 2010;Maliva et al. 2019;van Ginkel et al. 2014;Ward et al. 2009), the associated thinning of the injected freshwater volume has a negative impact on the RE (max 69% if the ambient fluid was brackish water versus around 89% if the ambient fluid was freshwater for the singlecycle experiments).
The vertical wells in the tank formed conduits for the freshwater that was injected in the bottom half of the aquifer. The buoyancy of the freshwater caused it to migrate through the vertical tubes, which had a much higher permeability than the surrounding glass beads. Such short-circuiting along wells has also been observed in operational MPPW-ASR systems Fig. 9 Change of the recovery efficiency (RE) with multiple cycles Fig. 10 Scatter plot showing the model-calculated versus the measured recovery efficiencies  and plays a role in karst aquifers where freshwater moves preferentially along natural vertical conduits (Missimer et al. 2002). The possibility of undesirable redistribution of freshwater through ASR well infrastructure by this process should thus be considered when investigating ASR construction options to optimize system performance in brackish aquifers such as separate wells with long screens (Maliva et al. 2006). During the recovery phase, the tubes could also act as a conduit for the upward moving brackish water.
Experiments in which freshwater occupied the entire thickness of the aquifer early on during the injection phase where not susceptible to the formation of salt fingers. This was the case for the vertical wells and the four point wells under pumping scheme 2, but for the two point wells and pumping scheme 1, gravitational instabilities were observed. Pumping scheme 2 is thus preferable as it prevents the formation of a disperse zone of intermediate salinities along the top of the freshwater, which may negatively impact the RE.
Up-coning below the wells is another cause for the relatively low RE in brackish aquifers. The difference between experiments 6 and 12 during the recovery phase exemplifies how the abstraction of brackish water from the deeper well prevented up-coning of the brackish water. This is the prime reason why the RE of pumping regime 2 was systematically greater than that of pumping regime 1, which reinforces findings from other studies that that scavenger wells can greatly enhance the performance of MPPW-ASR systems (Stuyfzand and Raat 2010;van Ginkel et al. 2014).
The injection and recovery rates also critically influenced the RE. The RE was consistently lower for experiments in category L than category H, due to the increased importance of the buoyancy effect, which is expected as forced convection (the flow caused by pumping) decreases relative to the free convection (the flow due to density differences) by a factor of 2 (Ward et al. 2007). The buoyancy results in a narrower shape during the injection phase and it accelerates the thinning of the freshwater volume during the recovery phase (Fig.  7). In a field-operational setting this means that minimizing the storage cycle duration as well as maximizing the abstraction rate would enhance the performance, but since the timing of the demand controls how much water needs to be provided and when, this is a system variable that is not easily controlled.
In previous studies, the RE of brackish ASR sites has been shown to improve with increasing number of injection and recovery cycles (Merritt 1986;Pyne 2015;) as the unrecovered water from the previous cycle forms a buffer zone along the perimeter of the injected freshwater. In the multicycle experiments of this study, the formation of such a buffer zone was clearly visible, albeit that it occupied the top part of the aquifer only-see time-lapse videos: (ESM2-ESM8). This confirms the assertion by Zuurbier et al. (2014) that in the deeper part of the system, the starting conditions for later cycles remain very similar to those during the first cycle. The presence of brackish water at a short distance below the abstraction well means that the system is vulnerable to salinization by up-coning. Nevertheless, the RE increased with every cycle (Fig. 9).

Study limitations
The setup of the laboratory tank forms a limitation in that the 3D radially dominated flow patterns, as expected under field conditions, are not simulated. Nevertheless, this type of experiment forms an important addition to the suite of available investigation techniques, which further includes field trials or mathematical modelling. The laboratory experiments have a big advantage over field trials, in which concentration patterns can only be observed indirectly (Haaken et al. 2017) and which are inherently associated with incomplete knowledge of aquifer properties (Miotliński et al. 2014;Zuurbier et al. 2014).
In experiments where brackish water was sitting on top of freshwater, salt fingers were observed in the laboratory tank (Fig. 6). These did not form in the numerical model, but instead a rather irregular salinity pattern was simulated by the numerical model, which did roughly coincide with the outline of the region of salt fingering as observed in the tank. The failure of the model to produce salt fingers is tentatively attributed to the fact that dispersive mixing in the simulations was overestimated, despite the very low dispersivity used, thus suppressing the instabilities that lead to the formation of salt fingers.
The experiments overestimate the RE that are reported for MPPW-ASR systems under field conditions (Zuurbier et al. 2014). One reason is that the effect of a background groundwater flow field, which causes the injected water to drift away from the well (Ward et al. 2009;, was not considered. Moreover, freshwater may be lost by upward seepage through aquitards (Zuurbier et al. 2014), whereas the model aquifer considered in this study was fully confined. Another important reason is that the glass beads form a porous medium that is much more homogeneous than real-world aquifers. In natural target aquifers, heterogeneities cause mixing between the injected freshwater and the ambient freshwater, resulting in a large water volume along the fringes in which salinities exceed the water quality threshold values (Miotliński et al. 2014). On the other hand, the laboratory experiments are likely to overestimate the negative impact of the freshwater buoyancy on the RE, because of the low anisotropy of the beads.

Conclusions
This study consisted of a series of laboratory tank and corresponding numerical modelling experiments to investigate the fate of freshwater injected into, and abstracted from, brackish aquifers using MPPW-ASR technology. The aim of the study was to determine how well design, pumping rate and pumping scheme influence the RE of MPPW-ASR systems in a homogenous, confined aquifer. This was achieved by implementing 16 combinations of different pumping rates (high or low), pumping schemes (with or without scavenger wells) and well designs (vertical wells with long screens, 2 or 4 partially penetrating wells) in single cycle experiments, followed by three multicycle experiments testing the development of the RE in three set-ups.
The visualization by time-lapse photography is the only means to directly observe the redistribution of the freshwater during the injection, storage and recovery cycles. Limitations of the laboratory conditions compared to the field are the lack of 3D flow, the absence of interfering factors that influence the RE like background groundwater flow and heterogeneities, and a lower vertical anisotropy. The calculated REs are therefore higher than what may be expected for real-world systems. Nonetheless, the insights gained highlight matters of concern in the planning phase and assist in developing conceptual understanding based on the interpretation of field data. In particular, the experiments conducted in this study allow making the following observations about the effect of MPPW-ASR system design on the RE: & Lower injection and abstraction rates (category L experiments) increased the duration of the injection and recovery cycles, providing a longer time window for upward drift of the freshwater, which reduced the RE compared to high injection/abstraction rates (category H experiments). This finding reinforces the importance of short cycle durations for successful ASR in brackish aquifers. As cycle lengths in in-situ systems usually are determined by the phases of excess water availability and water demand, the possibilities of making adjustments in existing systems are limited, but it must be a prominent consideration in the planning process, including MAR suitability mapping and decision support tools. & For the category L experiments, there were no significant differences in performance between the well systems considered. For the experiments in category H, however, the four-point well system mostly performed better than the vertical wells or the two point wells. The advantage of a larger number of screens in different depths seems to depend on the cycle duration. As the construction and operating costs of in-situ systems increase with the number of screens/wells, modelling different well configurations prior to implementation could help reduce unnecessary costs.
& The performance of the two-point well system was particularly poor when (1) an unstable density stratification (i.e. brackish water above freshwater) developed during injection, resulting in enhanced mixing by salt fingering, and (2) when no scavenger well was used. In in-situ systems, the first issue can be addressed by injecting a proportion of the freshwater through the shallow screens or wells. The application of scavenger wells is a more complex issue of economic viability, as the treatment costs or costs of discharge have to be balanced against the potential benefits. & Even though the RE of the multicycle experiments increased with successive cycles, the MPPW-ASR system remained susceptible to salinization. This is because the freshwater buffer zone that forms is essentially just a thin veneer along the aquifer top. Without a freshwater buffer zone below the abstraction well, there is the risk of salinization by saltwater up-coning. Well design and pump operation that minimise drawdown in the uppermost part of the aquifer, and thus the propensity for upward flow, are therefore advisable. & Freshwater could rapidly flow to the top of the aquifer during injection through conduits formed by the vertical wells. During recovery the vertical wells formed preferential pathways for upward flow of brackish water. The presence of multiple wells can thus form a risk for a scheme's performance. For in-situ systems, multi-level wells (with multiple screens in one borehole) could reduce the risk of preferential flow, provided hydraulic short-circuiting through the borehole annulus can be prevented. This is technically challenging and might not be feasible everywhere.
Finally, it is noted that the numerical model's ability to reproduce the visible salinity distribution in the tank, as well as the measured REs, confirms that they are suitable to study the performance of MPPW-ASR systems. The challenge with model applications to field situations, however, still lies in the knowledge of local site conditions, which is always uncertain and incomplete. Nevertheless, numerical modelling, preferably in combination with more complex laboratory experiments or field trials, provides a way to study MPPW-ASR performance for an extrapolated range of conditions.