Abstract
This paper explores post-pumping seawater intrusion (PP-SWI), which is the phenomenon of seawater intruding further inland than the location of a well, after pumping has ceased. Despite numerous papers on the topic of seawater intrusion and pumping, this is the first time that PP-SWI has been described in the literature, to our knowledge. This paper describes a laboratory-scale investigation of the phenomenon and we demonstrate that PP-SWI can be reproduced within physical experiments. We also show, using numerical modelling, that PP-SWI is caused by disequilibrium in the flow field following the cessation of pumping. Specifically, in our simulations, the cone of depression persisted after the cessation of pumping (first moving inland and then retreating toward the coastal boundary) which caused a lag in the reestablishment of fresh water flow toward the coast, after pumping had stopped. It was during this period of flow-field disequilibrium that PP-SWI occurred. We expect systems with larger postextraction disequilibrium to be most susceptible to PP-SWI and recommend future research to improve understanding of the relationship between hydrogeological parameters, extraction rates, well location, and incidence of PP-SWI. In those systems where PP-SWI is most likely, quantitative analysis of groundwater extraction and SWI will need to employ transient approaches to ensure that SWI is not underestimated.
1. Introduction
Coastal aquifers are major sources of freshwater supply in many countries [1–7]. Changes in the hydrology of the coastal zone, through groundwater pumping for example, can cause landward movement of seawater, a process referred to as seawater intrusion (SWI) [8]. SWI causes loss of water security in coastal aquifers through the displacement of fresh groundwater by seawater and through the degradation of groundwater quality whereby less than 1% of seawater (~250 mg/L chloride) renders freshwater unfit for drinking [9, 10].
SWI caused by groundwater pumping has received considerable attention over the years. In a recent review of the literature on field-based descriptions of well salinization and SWI, Houben and Post [11] noted the first systematic study to be that of Stephenson [12] who reported an increase in salt concentrations in water supply wells in Liverpool, England. This study predates by 50 years the famous paper of Herzberg [13], who described increasing salt concentrations in groundwater wells on the German tourist island of Norderney, following an unusually dry winter and spring. Herzberg [13] also described the relationship between freshwater and saltwater (i.e., the Ghyben-Herzberg approximation), which is fundamental to the quantitative modelling of SWI.
Despite the long history of studies detailing groundwater pumping and SWI, assessing the likelihood of groundwater withdrawals causing well salinization is still not straightforward. For example, the simulation of variable-density groundwater flow and solute transport in transient and heterogeneous coastal groundwater systems is computationally expensive, requires extensive field data for calibration, and is not always successful [14–16]. Data scarcity and faster computational times have meant that sharp interface analytic models are widely used to study SWI problems [e.g., 17, 18]. A well-known approach is that of Strack [19], who used potential theory to solve for the position of the freshwater-saltwater interface in regional groundwater flow systems with a pumping well. The solution also estimates the critical (i.e., maximum) pumping rate allowable before the well becomes salinized. The physics contained in analytical solutions and the simplicity of their application can provide a good balance between accounting for the relevant processes, data requirements, and calculation times [20, 21].
However, it is important to understand the implications of the simplifying assumptions underlying analytic solutions. One of the most important of these assumptions is that of steady-state conditions (i.e., the extent of seawater within the aquifer is estimated after infinite time). Research into SWI and sea-level rise (SLR) has shown that the steady-state assumption may lead to an underestimation of SWI impacts because of the so-called “seawater intrusion overshoot” phenomenon [22–26]. Seawater intrusion overshoot is the temporary movement of the saltwater-freshwater interface inland beyond the post-SLR steady-state position. Morgan et al. [25, 26] used physical and numerical modelling to show that the overshoot is physically reproducible and can occur within field-scale aquifers under realistic sea-level rise scenarios, as predicted by the Intergovernmental Panel on Climate Change [27]. Morgan et al. [26] found that the overshoot occurred during the period of flow-field disequilibrium (i.e., when the head and flow had not reached their post-SLR steady state) and was magnified when the disequilibrium was longer.
In this study, we consider a phenomenon that has not been described previously (to our knowledge) and that is likely also caused by flow-field disequilibrium, but in this case following the cessation of pumping. In a modelling study by Walther et al. [28] of the Oman coastal aquifer system (where groundwater has been intensively used for agricultural purposes for several decades under a semiarid climate conditions), scenario calculations showed that salinities in some of the wells did not immediately drop after pumping stopped. Rather, salinities in these wells increased for a certain amount of time (i.e., years to decades) suggesting that seawater continued to move landward beyond the pumping well location after the cessation of pumping. We term this phenomenon post-pumping seawater intrusion (PP-SWI) and hypothesize that disequilibrium in the flow field following pumping is causing the PP-SWI phenomenon.
In this study, we use a physical experiment in combination with numerical modelling to investigate the PP-SWI phenomenon in a systematic manner. Specifically, we aim to show that SWI does not necessarily stop immediately after the cessation of pumping in a coastal aquifer, and may penetrate further inland than the location of the pump after the cessation of pumping. Given this objective, laboratory conditions were designed to be especially conducive to PP-SWI. This included the use of a fixed-flux inland boundary condition, shown to increase flow-field disequilibrium in the seawater intrusion overshoot studies of Morgan et al. [25, 26]. Also, active seawater intrusion was instigated to enhance SWI within the study. This is because, under active SWI, the hydraulic gradient slopes from the ocean toward the well under pumping, as described previously by Badaruddin et al. [29] and Werner [30]. Unconfined conditions are considered because these types of aquifers are most commonly used for water supply purposes [22]. The present study neglects direct groundwater recharge from above, which is a valid assumption for arid regions, e.g., Arab countries and Australia.
2. Methods and Materials
2.1. Physical Model
An acrylic box of 2.0 m length, 0.5 m height, and 0.05 m width was used for physically simulating a cross section of an unconfined aquifer consisting of quartz sand (Figure 1). The apparatus used for physical modelling had essentially the same dimensions as that used by Stoeckl and Houben [31] and Strack et al. [18] in their investigation of freshwater lenses and saltwater intrusion. The box was homogeneously filled with quartz sand up to a height of 0.2 m. Using a Camsizer® (Retsch Technology, Germany), with a measurement range of 30 μm to 30 mm, the grain size distribution of the filter sand ( and ) was optically determined. Applying Hazen’s equation [32] on the grain size distribution, a hydraulic conductivity (L/T) with a value of was determined. A mean of was additionally obtained using Darcy’s constant-head conductivity tests. Resulting intrinsic permeabilities are and , respectively, considering a density of 996.9 kg m−3 and a viscosity of 0.00089 kg m-1 s-1, for freshwater at 25°C [32]. Porosity (-) was calculated by density measurements using a helium pycnometer (Micromeritics®, Germany) and was found to be 0.41. For visualization purposes, freshwater and saltwater were coloured using different tracer dyes, at a concentration of 0.3 g L−1. Freshwater was coloured with uranine (yellow) and indigotine (blue) to distinguish flow paths, whereas saltwater was coloured with eosin (red). Dyed freshwater and saltwater densities were determined using a DMA 38 density meter (Anton Paar, Austria) and were found to be 996.9 kg m−3 (freshwater) and 1020.9 kg m−3 (saltwater) at 25 °C.
Sixteen tubes were attached to a peristaltic pump, creating a fixed-flux inland boundary condition. Tube outlets were equally spaced in the vertical direction below the water table on the right-hand edge of the acrylic box. The flow was distributed equally between the tubes, with a total flux of 0.086 m3 d−1, which was maintained throughout the entire experiment. A constant water level was maintained in the saltwater reservoir at the left side of the sand box, simulating the sea. A permeable fabric (Huesker Synthetic, Germany) with a mesh spacing of 0.45 mm was used to construct a vertical boundary, representing the coastline. According to the physical properties of the fabric, i.e., mesh spacing, no additional resistance in water flow is generated. Freshwater discharge from the aquifer into the “ocean” formed a thin layer in the saltwater reservoir, which was removed by an overflowing bucket at a tube’s fixed position. The saltwater reservoir was continuously recharged with additional saltwater to maintain the water level and avoid dilution.
A tube with a 2 cm perforation at the base was placed 60 cm from the simulated coastline at a height of 4 cm above the aquifer base (i.e., the base of the sand tank), as shown in Figure 1. The tube was connected to another peristaltic pump such that water could be extracted and a well simulated in the physical model. Initially, constant boundary conditions were maintained for around 12 h to achieve a steady-state interface. Steady-state conditions were then confirmed by monitoring the interface position for a period of 1 h, over which time the interface remained stationary. The interface was visually approximated as the midpoint of the transition zone between freshwater and saltwater in the aquifer. Groundwater pumping was simulated by extracting water from the well at a constant rate of 1 m3 d−1. Water extraction was stopped just prior to the interface reaching the well, and the experiment continued until the interface returned to the original steady-state position. The entire experiment was filmed at an interval of one picture per second (Sony XDCAM EX) enabling detailed optical analyses of quantitative indicators described in Section 2.3.
2.2. Numerical Modelling
The physical experiment was simulated using the finite-element code FEFLOW 7.0, which considers variably saturated, density-dependent flow and solute transport processes [33] and is benchmarked for variable-density flow by Stoeckl et al. [34]. The numerical model was assumed to be a unit-width aquifer (2-D), and hence, the freshwater inflow at the inland boundary was scaled to .
The vertical coast was simulated using the Dirichlet boundary conditions for head and concentration (Figure 2). A constraint was applied to the Dirichlet-type concentration boundary such that seawater concentration () occurred only where inflow occurs. The concentration of inflows through the inland boundary on the right-hand side was (representing freshwater). Inflows at the inland boundary were applied using a Neumann boundary condition. As in Morgan et al. [25], a seepage face boundary condition was implemented in the numerical model along the vertical beach face. The seepage face boundary was implemented using a fixed-head boundary (with the head equal to the elevation of each node) that is only active if there is flow out of the model. Otherwise, the vertical beach face was a no-flow boundary above the seepage face. A mass sink term was implemented at one single node, i.e., 0.6 m inland from the coastal boundary and 0.05 m above ground, representing the center of the well screening in the physical model.
A combined quadrilateral-trapezoidal mesh with around 94,500 elements was used. The model grid was refined in the area between the coastal boundary and the well to allow for more stable calculations. Grid discretization was , and for the refined area 1 mm. Automatic time stepping was used with an upper step size limit of 0.02 s.
The unsaturated zone was modelled using van Genuchten’s [35] functions, with different curve fitting parameters for alpha (L-1) and (-). Initial values for and were taken from the physical experiment and varied during calibration, together with van Genuchten’s [35] parameters. A density ratio (-) of 0.0241 was applied in accordance with the physical modelling values, defined as , where and are freshwater and saltwater densities (mL−3), respectively. Longitudinal and transverse dispersivities were set to and , respectively, and the molecular diffusion coefficient was set to .
2.3. Quantitative Indicators
Several indicators were used to measure saltwater intrusion and retreat, as well as water level in the vicinity of the well. The wedge toe occurs at the point of intersection between the aquifer base and the freshwater-saltwater interface. The wedge toe length is defined as the distance between the vertical coastal boundary and the wedge toe. In the physical experiments, the interface was visually determined as the transition between the freshwater (yellow/blue in colour) and the saltwater (red in colour), as described above. In the numerical modelling, the interface was characterized as the 50% saltwater concentration isochlor.
The cone of depression was also used as a quantitative indicator. The cone of depression was visible as a grey-shaded area above the pump in the physical model, when pumping was enabled. For visualization purposes and for a comparison to numerical model results, a blue line was drawn onto the photographs by hand, indicating the location and extent of the cone of depression more clearly. In the numerical model, the cone of depression was visualized using a grey scale for moisture content between 0.25 (unsaturated) and 1 (fully saturated).
Freshwater flow paths were also used as indicators and were visualized by two different tracer dyes in the physical model. The changes in the flow paths were tracked over time and provided additional information relating to the transient response to pumping within the freshwater component of the aquifer. In the numerical model, velocity contours and vectors were used as indicators of flow direction for the freshwater and saltwater.
3. Results and Discussion
3.1. Physical Model
A steady-state equilibrium between saltwater and freshwater was initially established in the two-dimensional physical experiment, without pumping. At steady-state, the saltwater wedge toe length was 0.09 m due to the density contrast between freshwater and saltwater. This equilibrium condition can be seen in Figure 3(a) (). Under equilibrium conditions, flow paths can be seen to be roughly parallel to the aquifer base, except above the saltwater wedge where they are bent upwards, as expected.
The pump was switched on at time , and the saltwater wedge started to move inland (compare Figure 3(b), and Figure 3(c), ). The pump was switched off at time , just before seawater entered the well. A valve stopped the water in the well from flowing back into the sand tank. After pumping had stopped, the saltwater continued moving inland until its wedge toe length reached a maximum value of 0.75 m at . That is, after the cessation of pumping, the saline interface continued to move inland beyond the pump location and at was 0.15 m further landward than the pump location. After the maximum wedge toe length was reached at , the wedge toe slowly started receding back to its original steady-state position, which was reached after around 5000 s.
Immediately after the pumping commenced (at ), the yellow and blue freshwater flow paths bent toward the pumping well (see Figure 3(b), and Figure 3(c), ), as expected. After pumping stopped (at ), this zone of disturbed flow paths can be seen to move inland beyond the pumping well, indicating the inland movement of freshwater and saltwater. That is, freshwater as well as saltwater (within the saltwater wedge) was observed to move landward following the cessation of pumping. After around , the zone of disturbed freshwater flow paths was observed to start retreating back toward the coast.
A cone of depression was identified as a result of pumping, as expected. The position and extent were tracked over time in the video that was recorded during the experiment. In Figure 3, the position of the saturated-unsaturated interface was drawn in by hand for each of the selected time frames. The cone of depression started forming immediately after the pump was switched on and continuously extended vertically as well as horizontally until it almost reached the depth of the pump location at , at which time the pump was turned off. The pump did not take in air at any point in time, which confirms that the maximum depth of the cone of depression always remained above pump level. When the pump was switched off, the cone of depression immediately started to recover. As described previously for the zone of freshwater flow disturbance, the cone of depression moved landward after the cessation of pumping before retreating back toward the coast.
3.2. Numerical Model
The establishment of a stable simulation of the “small scale” sand tank experiment incorporating density driven flow, as well as unsaturated flow, was time consuming and difficult to achieve. The mesh was refined several times, especially in the region between the coast and the pump. The single pump location was eventually also redistributed to 10 different “well nodes,” each abstracting 1/10th of the total volume. Time steps had to be as low as 0.01 s to enable stability. In addition to the testing of grid resolution and time-stepping schemes, different solvers (e.g., pcg, pardiso, and samg) were tested and the samg solver showed the best performance and stability; it was subsequently used for the simulations. The calibration was done manually, because the simulations had to be stopped and restarted to enable the switching on and off of the groundwater well for extraction (compare [22]). Every simulation was done in three steps: (a) achieving an initial steady state, (b) enabling groundwater extraction by implementing a well, and (c) disabling groundwater extraction by removing the well. The first and third phase took 5000 s each, whereas the second phase (where groundwater is being extracted) was set to 135 s, in line with the period of time that groundwater was extracted in the physical experiment.
During calibration of the numerical model, we tried to match several observations from the physical experiment: (a)Steady-state wedge toe length of 0.09 m at (b)Wedge toe length of 0.6 m (i.e., at the pump location) at (c)Maximum wedge toe length of 0.75 m at (d)A similar time (of 5000 s) for the wedge toe to return to the preextraction steady-state location(e)A realistic saltwater wedge shape at all times that was comparable to the observations in the physical experiment
Parameters taken into account during the calibration process were (varied between and ), (between 0.25 and 0.4), and van Genuchten’s [35] parameters alpha (varied between 1.3 m-1 and 20 m-1) and (varied between 0.3 and 3). As the first two parameters were measured in the laboratory, they were varied over a smaller range than the latter two. It was observed that the variation of van Genuchten’s [35] parameters had a strong influence on the temporal as well as the spatial evolution of the saltwater wedge.
Out of the approximately 150 different parameter combinations trialled, simulated wedge toe trends from a selection of 15 parameter combinations are shown in Figure 4. The toe length trend in all of these simulations reflect the typical behaviour of saltwater intruding due to pumping followed by recession due to the cessation of pumping, as observed in the physical experiment (black dots in Figure 4). However, many parameter sets do not show a temporary continuation of saltwater intrusion after the cessation of pumping, some do not show an overall intrusion length further than the pump location at 0.6 m, and for some parameter combinations the initial and final steady-state wedge toe positions differ from those observed in the physical model (of 0.09 m). One parameter set with a porosity of 0.3, van Genuchten’s [35] parameters and , and of was selected for detailed analysis (black line in Figure 4). This parameter set satisfied all but one of the previously mentioned criteria (a) to (e). The time for the saltwater wedge to return to its preextraction steady-state location was better met by using a porosity of 0.4 (compare a similar curve to the selected fit in Figure 4). However, this scenario deviates more from the observed maximum wedge toe length of 0.75 m at around . More importantly, a higher porosity leads to a quicker advance of the saltwater-freshwater interface during pumping and therefore to a significant intake of saltwater into the pump. This results in an awkward deformation of the saltwater wedge (not shown). Thus, the parameter set with the lower porosity of 0.3 was chosen to work with in the following model, even though the time of recession did not exactly match the observations in the physical experiment. The shorter recession times in the numerical model may have different causes: (a) hysteresis effects in the unsaturated zone may slow down the reequilibration process, as the drying process is faster than rewetting (the van Genuchten relationship of water suction capacity and water content in the numerical simulation does not account for hysteresis); (b) a significant portion of the sand tank is unsaturated at pump stop, thus certain pores may not have fully resaturated (due to entrapped air), which possibly leads to a reduced permeability and therefore retention in the backward movement of the saltwater-freshwater interface; and (c) the distribution of smaller and larger pores (double porosity) may additionally lead to a longer tailing (and blurring), as the coloured sea water is slowly released from the smaller pores. As the numerical model does not account for these effects, shorter times of recession are likely, compared to the physical experiment or even to real-world observations. For the interpretation of the results, however, the exact timing of the recession is less important than accurately simulating the interface movement to, beyond, and back to the pump location at 0.6 m, which is defined as PP-SWI here.
Interface location results for the selected fit simulation are shown in Figure 5. In addition to the 0.5 relative salinity contour, we also show the 0.1, 0.3, 0.7, and 0.9 relative saltwater contours as they provide information on the interface width. Between the commencement of pumping (at ) and the cessation of pumping (at ), the interface moved inland, as expected. At exactly , the 0.1 relative salinity contour reaches the well location, which was also the time when the pump was switched off in the physical experiment. The cone of depression (i.e., the white area above the pump) at this point is the largest. A widening of the interface was also observed and continued until a maximum wedge toe length of 0.76 cm was simulated at (compared to the physically measured time of ). The wedge toe receded back to the initial steady-state position faster in the numerical modelling than was observed in the physical experiment, as noted above.
During the wedge toe recession phase of the simulation, the wedge toe was temporarily closer to the coast than the initial steady-state position. After the wedge toe reached its minimum length of around 0.07 m, it then proceeded inland to reach its final steady-state position at 0.09 m. This overshooting of the steady-state position is visible in the last two time frames in Figure 5, where it can be seen that the 0.5 relative saltwater contour is closer to the coast at than at . The overshoot is also evident in Figure 4, where the line of the “best fit” simulation drops below the prepumping steady-state wedge toe length of 0.09 m (indicated by a dashed line) and then returns to the steady-state position at around . The wedge toe recession overshoot was not detected in the physical experiment, perhaps due to a blurring of the tracer dye and the small scale of the effect. The overshoot during the wedge toe recession (in response to a sea-level drop) was observed previously by Morgan et al. [25] using both physical and numerical modelling. As mentioned previously, the cause of the overshoot was shown by Morgan et al. [25, 26] to be disequilibrium between the hydraulic head and the saltwater interface as they both moved to a new steady state in response to a change in sea level. The overshoot observed in the current numerical modelling can also be attributed to this form of disequilibrium.
In Figure 5, grey shading is used to indicate the relative degree of saturation. Light grey-shaded areas at the top of the model indicate areas of lower saturation while black-shaded areas are saturated. It can be seen that at , there is a gradient in water saturation from right to left, which conforms to the flow of freshwater from right to left within the model. At (i.e., 10 s after pumping has commenced) a cone of depression had already developed above the well location. The cone of depression is visible as an area of lower relative saturation and is slightly asymmetric, with lower saturation over a greater depth inland of the well, compared to the coastal side of the well. The cone of depression gets larger until the cessation of pumping at , after which time the cone of depression slowly recovers. The relative saturation contours (which also indicate the approximate water table location and hence hydraulic gradient, not accounting for density effects) show that the cone of depression is causing saltwater and freshwater to be pulled toward the well from the coastal and landward boundaries, respectively. A link between the cone of depression and the direction of the saltwater movement (i.e., intrusion versus retreat) was observed. After the cone of depression was “filled up” by both the intruding seawater from the coast boundary and the freshwater from the landward boundary, the wedge toe started to retreat (at ). That is, the wedge toe began retreating once the cone of depression had recovered enough such that the hydraulic gradient and flow toward the coast in the vicinity of the saltwater wedge had reestablished, at . At , the steady-state relative saturation condition was reestablished with a saturation gradient from the land toward the ocean.
In Figure 6, Darcy flow velocities are visualized as black contours of 0.1, 0.5, 1, 2, 5, 10, and 20 m/s. At steady state, flow directions are generally horizontal from right to left. Freshwater is only diverted upward above the saltwater wedge circulation cell, as previously observed in various studies [34, 36, 37]. As pumping commences, the saltwater wedge starts intruding inland and flow velocity vectors are directed toward the pump (not shown here). Flow velocity contours in general stay parallel until pumping is stopped at (compare Figure 5). After , semicircular flow velocity contours start to develop at the bottom of the model, to the right of the well location. These are contours of low flow velocities and they indicate a singularity, where flow from the right and left sides merge and get directed upward. An excerpt view of the singularity (from Figure 6, time frame ) is shown in Figure 7 with additional flow vectors as black arrows. The full height of the numerical model is displayed, a scale in the upper-left corner shows the exact dimensions. Flow at the bottom of the model merges at the singularity from both sides where velocities are lowest (semicircular contours), and is directed vertically upward to fill the cone of depression above the singularity.
The singularity (and cone of depression above it) moved inland first (Figure 6, to ), before receding back toward the pump location (Figure 6, to ). As the semicircular contours surrounding the singularity reach the position of the intruding saltwater wedge at , they dissolve and form parallel lines again, until the initial condition is reached for both the concentration and Darcy velocity contours.
The singularity (and cone of depression) moved in and out of the aquifer faster than the saltwater wedge (Figure 6). We hypothesize that the occurrence and duration of this singularity dictate the extent to which seawater will continue to intrude after the cessation of pumping. Future studies are needed to investigate this further and to explore the influence of hydraulic parameters on the dynamics of the singularity. We expect that it will not be easy to draw direct linear relationships between single parameters and the extent of PP-SWI. It is rather expected that a local maximum for each parameter (or combination of parameters) exists, where the occurrence of the PP-SWI phenomenon is greatest. For example, when hydraulic conductivity is low, seawater intrusion will be slower and hence the likelihood of the PP-SWI phenomenon will be lower; however, if hydraulic conductivity is high, the cone of depression will develop more slowly and this might also reduce the likelihood of the PP-SWI phenomenon. We also expect that pump distance from the coast and extraction rate will influence the extent of PP-SWI in a nonlinear manner.
The upscaling of our findings to field-scale cases is not straight forward, as inherent difficulties in scaling sand box models exist [38]. Further discrepancies in transferring our findings to real-world cases result from neglecting, e.g., heterogeneity in the subsurface, direct groundwater recharge from above, or variable pump schemes in the sand tank experiment. The finding that PP-SWI is physically based, however, may be transferred from the two-dimensional setting used in the present study into the third dimension by considering a coastline, where pumping wells (e.g., well fields) are located a certain distance along the coastline and groundwater flow can be considered perpendicular to the coast. It is assumed that PP-SWI occurs at a field scale on much larger time scales. Here too, larger areas will be salinized that would be the case without PP-SWI, potentially threatening inland ecosystems and fresh groundwater resources.
PP-SWI might also be expected in certain cases where countermeasures are taken, for example, where pump rates are reduced (but not ceased), as well as when seawater has not reached the well but is still moving toward it. This is expected for cases where pump rates are not sustainable and groundwater tables are strongly lowered due to overexploitation over several decades, resulting in a large regional cone of depression in close vicinity to the coast. Times for filling up such a deficit water volume depend on recharge rates, as well as aquifer geometries. It can be stated, though, that the longer the time for reestablishing water tables and gradients from the land toward the sea, the longer SWI will take place and the more pronounced PP-SWI is expected to be. A proper analysis of field-scale cases is beyond the scope of this paper, which is about the mechanistic understanding of the PP-SWI phenomenon that was previously observed in the numerical model simulation of a coastal aquifer in Oman [28].
4. Conclusions
The physical experiment conducted here highlights the PP-SWI phenomenon, whereby seawater intrudes further inland than the location of a well, after pumping has ceased. A numerical model with the same extent as the physical model was calibrated and showed a very similar behaviour with minor differences in the overall time duration of the phenomenon. The calibration was considered adequate with a parameter set which satisfied four out of five predefined criteria: the steady-state wedge toe length without pumping, the time for the wedge toe to reach the well location after pumping had commenced, the maximum wedge toe length, and the shape of the wedge. The only predefined criteria that was poorly met was the timing of saltwater retreat.
The modelling carried out within this study offers insight into the cause of the PP-SWI phenomenon, at least for the conditions used in this study (i.e., active SWI and a fixed-flux inland boundary). In both the physical and numerical models, the driving force for the PP-SWI phenomenon was identified to be the presence, duration, and movement of the cone of depression. As a result of pumping, water entered the model from both sides, from the coastal boundary as seawater and from the inland boundary as freshwater, to fill up the “deficit volume” of the cone of depression. Under pumping, the new steady-state position of the saltwater wedge occurred at the pump and did not intrude further inland because the pumping well is the point where all waters merge. When the pump was turned off, the cone of depression moved inland (as did the singularity beneath the cone of depression), which caused water to flow from the coast and inland in order to fill the cone of depression. As a result, the saltwater wedge continued to move landward after the cessation of pumping until the cone of depression had recovered enough that a hydraulic gradient (and fresh water flow) toward the coast had reestablished. As such, it has been demonstrated that the PP-SWI phenomenon is caused by pumping-induced transience in the hydraulic gradient between the land and the ocean. The transience of the hydraulic gradient is linked to the dynamics of the cone of depression (and the associated flow singularity) following the cessation of pumping.
The PP-SWI phenomenon has not, to the best of our knowledge, been mentioned previously in the literature, except for the Oman case [28]. The observation of seawater intruding further inland than the location of an abandoned well field led us to investigate the possibility of the physical occurrence of PP-SWI and its driving mechanism. In this study, we have shown PP-SWI to be a physically reproducible phenomenon. It has important implications for coastal groundwater aquifers all around the world and shows that seawater intrusion induced by pumping (a) may not immediately stop or recede when pumping is ceased, (b) may significantly extend further inland than the actual position of a previously active pump or well field, and (c) may have a direct impact on the time scales for remediating salinized wells. This means that not only will there be salinization of much greater areas of coastal aquifers, but the times for sea water to retreat will also be much longer than expected. In (semi-)arid areas with limited groundwater recharge, such as in Arab countries, in Australia, or in other coastal areas of Africa and South America, SWI associated with groundwater pumping may have been underestimated if transient processes were not accounted for. SWI induced by freshwater abstraction from coastal aquifers is not limited to the time period of pumping, but rather extends beyond it, depending on pumping rates and locations, and aquifer characteristics, e.g., hydraulic conductivity and porosity, length-to-height relationship of an aquifer, density contrast between freshwater and saltwater, and recharge condition. It is thus important to foster research of transient dynamics in groundwater systems rather than to assume steady-state solutions, especially when investigating the effects of groundwater abstraction in coastal aquifers.
Data Availability
Measurement data of the saltwater-freshwater interface were generated by the experimental work and are documented in the manuscript.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors thank Katrin Damm for helping with the physical experiments. Leanne Morgan is partly supported by Environment Canterbury, New Zealand (https://www.ecan.govt.nz/).