Aquifer Storage and Recovery (ASR) is a managed aquifer recharge method where water is injected and later extracted using wells. In saline aquifers, ASR performance is often limited by dispersive mixing, which creates a transition zone at the edge of the injected freshwater and buoyancy-driven flow, which causes the freshwater to rise and deform during storage—both reducing recovery efficiency. This study investigates whether horizontal wells can improve ASR performance in saline, low-transmissivity aquifers by achieving acceptable recovery efficiencies and outperforming conventional vertical wells. Three configurations were evaluated numerically with MODFLOW 6: a horizontal well, a fully penetrating vertical well, and a dual well system with a fully penetrating injection well and a partially penetrating extraction well. Models were tested on a large set of parameter combinations from Latin Hypercube Sampling, targeting conditions where vertical wells perform poorly. The horizontal well generally achieved higher recovery efficiencies, with a median of 45% after five ASR cycles, compared to 6% and 16% for the fully and partially penetrating vertical wells. Its advantage was greatest under strong buoyancy conditions, where vertical wells failed to recover any freshwater. While dispersive mixing reduced horizontal well performance by causing earlier saltwater breakthrough, it improved vertical well recovery by stabilizing the injected freshwater. In conclusion, horizontal wells are promising for ASR when hydraulic conditions require multiple vertical wells and when buoyancy-driven flow significantly limits vertical well performance.

Lin and Lin (Faulty assumptions: Groundwater modeling through anisotropic fault zones, Journal of Hydrology 653(2025)) make unfounded claims about earlier studies of fault aquifer interaction: they state that the standard boundary conditions for a conductive fault used in earlier studies are flawed or based on faulty assumptions and that the earlier studies did not consider the discontinuities in both the normal component of flow and hydraulic head across a general fault. Further, they are unclear about the approximations in their own analysis which attempts to replace the two-dimensional flow field within a thin, anisotropic fault zone with two, one-dimensional, internal boundary conditions. Their claims about earlier studies are refuted and the approximations in their analysis are examined. Despite the limitations of their analysis, the work of Lin and Lin has value.

Research into land subsidence caused by groundwater withdrawal is hindered by the availability of measured heads, subsidence, and forcings. In this paper, a parsimonious, linked data-driven and physics-based approach is introduced to simulate pumping-induced subsidence; the approach is intended to be applied at observation well nests. Time series analysis using response functions is applied to simulate heads in aquifers. The heads in the clay layers are simulated with a one-dimensional diffusion model, using the heads in the aquifers as boundary conditions. Finally, simulated heads in the layers are used to model land subsidence. The developed approach is applied to the city of Bangkok, Thailand, where relatively short time series of head and subsidence measurements are available at or near 23 well nests; an estimate of basin-wide pumping is available for a longer period. Despite the data scarcity, data-driven time series models at observation wells successfully simulate groundwater dynamics in aquifers with an average root mean square error (RMSE) of 2.8 m, relative to an average total range of 21 m. Simulated subsidence matches sparse (and sometimes very noisy) land subsidence measurements reasonably well with an average RMSE of 1.6 cm/year, relative to an average total range of 5.4 cm/year. Performance is not good at eight out of 23 locations, most likely because basin-wide pumping is not representative of localized pumping. Overall, this study demonstrates the potential of a parsimonious, linked data-driven, and physics-based approach to model pumping-induced subsidence in areas with limited data.

This paper presents the results of the 2022 Groundwater Time Series Modelling Challenge, where 15 teams from different institutes applied various data-driven models to simulate hydraulic-head time series at four monitoring wells. Three of the wells were located in Europe and one was located in the USA in different hydrogeological settings in temperate, continental, or subarctic climates. Participants were provided with approximately 15 years of measured heads at (almost) regular time intervals and daily measurements of weather data starting some 10 years prior to the first head measurements and extending around 5 years after the last head measurement. The participants were asked to simulate the measured heads (the calibration period), to provide a prediction for around 5 years after the last measurement (the validation period for which weather data were provided but not head measurements), and to include an uncertainty estimate. Three different groups of models were identified among the submissions: lumped-parameter models (three teams), machine learning models (four teams), and deep learning models (eight teams). Lumped-parameter models apply relatively simple response functions with few parameters, while the artificial intelligence models used models of varying complexity, generally with more parameters and more input, including input engineered from the provided data (e.g. multi-day averages). The models were evaluated on their performance in simulating the heads in the calibration period and in predicting the heads in the validation period. Different metrics were used to assess performance, including metrics for average relative fit, average absolute fit, fit of extreme (high or low) heads, and the coverage of the uncertainty interval. For all wells, reasonable performance was obtained by at least one team from each of the three groups. However, the performance was not consistent across submissions within each group, which implies that the application of each method to individual sites requires significant effort and experience. In particular, estimates of the uncertainty interval varied widely between teams, although some teams submitted confidence intervals rather than prediction intervals. There was not one team, let alone one method, that performed best for all wells and all performance metrics. Four of the main takeaways from the model comparison are as follows: (1) lumped-parameter models generally performed as well as artificial intelligence models, which means they capture the fundamental behaviour of the system with only a few parameters. (2) Artificial intelligence models were able to simulate extremes beyond the observed conditions, which is contrary to some persistent beliefs about these methods. (3) No overfitting was observed in any of the models, including in the models with many parameters, as performance in the validation period was generally only a bit lower than in the calibration period, which is evidence of appropriate application of the different models. (4) The presented simulations are the combined results of the applied method and the choices made by the modeller(s), which was especially visible in the performance range of the deep learning methods; underperformance does not necessarily reflect deficiencies of any of the models. In conclusion, the challenge was a successful initiative to compare different models and learn from each other. Future challenges are needed to investigate, for example, the performance of models in more variable climatic settings to simulate head series with significant gaps or to estimate the effect of drought periods.

The performance of time series models is assessed using synthetic head series simulated with a numerical model that solves Richards' equation for variably saturated flow. Heads were simulated in a homogeneous unconfined aquifer between two parallel canals; measured daily precipitation and potential evaporation are specified at the land surface and root water uptake is simulated. The head response to a precipitation event is nonlinear and depends on the saturation degree and rainfall before and after the precipitation event while evaporation reduction occurs during summers. Synthetic series were generated for 27 years and three different soil types; the unsaturated zone thickness varies between 0 and >5 m. The synthetic head series were simulated with a linear and nonlinear time series model. Performance of a linear time series model with four parameters, using a scaled Gamma response, gave R² values ranging from 0.67 to 0.96. The nonlinear time series model with five parameters simulates recharge using a root zone reservoir after which the head response to recharge is simulated with a scaled Gamma response function. The nonlinear time series model was able to simulate all synthetic head series very well with R² values above 0.9 for almost all models. The head response of the nonlinear model to a single precipitation event compares well to the response of the variably saturated groundwater model. The provided scripts may be used to simulate synthetic head series for other climates or for systems with additional complexity to assess the performance of other data-driven models.

Many sedimentary aquifers consist of small layers of coarser and finer material. When groundwater flow in these aquifers is modeled, the hydraulic conductivity may be simulated as homogeneous but anisotropic throughout the aquifer. In practice, the anisotropy factor, the ratio of the horizontal divided by the vertical hydraulic conductivity, is often set to 10. Here, numerical experiments are conducted to determine the effective anisotropy of an aquifer consisting of 400 horizontal layers of which the homogeneous and isotropic hydraulic conductivity varies over two orders of magnitude. Groundwater flow is simulated to a partially penetrating canal and a partially penetrating well. Numerical experiments are conducted for 1000 random realizations of the 400 layers, by varying the sequence of the layers, not their conductivity. It is demonstrated that the effective anisotropy of the homogeneous model is a model parameter that depends on the flow field. For example, the effective anisotropy for flow to a partially penetrating canal differs from the effective anisotropy for flow to a partially penetrating well in an aquifer consisting of the exact same 400 layers. The effective anisotropy also depends on the sequence of the layers. The effective anisotropy values of the 1000 realizations range from roughly 5 to 50 for the considered situations. A factor of 10 represents a median value (a reasonable value to start model calibration for the conductivity variations considered here). The median is similar to the equivalent anisotropy, defined as the arithmetic mean of the hydraulic conductivities divided by the harmonic mean.

Deformations of the earth’s crust create tortuous paths for groundwater flow, altering pressure distributions and flow lines. A solution for steady groundwater flow through a deformed aquifer is derived by applying the singular point method using a rectangular reference plane. The singular point method is used to develop conformal mappings with complex geometries using basic principles of groundwater mechanics including superposition, the method of images, and stagnation point analysis. The derived solution contains three parameters that can be chosen to simulate flow through a variety of deformed aquifers, including flow through a normal fault with a 90^∘ dip, flow through a fold, and flow through a relay ramp.

Global mean sea-level rise (SLR) has accelerated since 1900 from less than 2 mm yr-1 during most of the century to more than 3 mm yr-1 since 1993. Decision-makers in coastal countries, however, require information on SLR at the regional scale, where detection of an acceleration in SLR is difficult, because the long-term sea-level signal is obscured by large inter-annual variations with multi-year trends that are easily one order of magnitude larger than global mean values. Here, we developed a time series approach to determine whether regional SLR is accelerating based on tide gauge data. We applied the approach to eight 100-year records in the southern North Sea and detected, for the first time, a common breakpoint in the early 1990s. The mean SLR rate at the eight stations increases from 1.7 ± 0.3 mm yr-1 before the breakpoint to 2.7 ± 0.4 mm yr-1 after the breakpoint (95% confidence interval), which is unprecedented in the regional instrumental record. These findings are robust provided that the record starts before 1970 and ends after 2015. Our method may be applied to any coastal region with tidal records spanning at least 40 years, which means that vulnerable coastal communities still have time to accumulate the required time series as a basis for adaptation decisions in the second half of this century.

Slow-moving landslides move downslope at velocities that range from mm year⁻¹ to m year⁻¹. Such deformations can be measured using satellite-based synthetic aperture radar interferometry (InSAR). We developed a new method to systematically detect and quantify accelerations and decelerations of slowly deforming areas using InSAR displacement time series. The displacement time series are filtered using an outlier detector and subsequently piecewise linear functions are fitted to identify changes in the displacement rate (i.e., accelerations or decelerations). Grouped accelerations and decelerations are inventoried as indicators of potential unstable areas. We tested and refined our new method using a high-quality dataset from the Mud Creek landslide, CA, USA. Our method detects accelerations and decelerations that coincide with those previously detected by manual examination. Second, we tested our method in the region around the Mazar dam and reservoir in Southeast Ecuador, where the time series data were of considerably lower quality. We detected accelerations and decelerations occurring during the entire study period near and upslope of the reservoir. Application of our method results in a wealth of information on the dynamics of the surface displacement of hillslopes and provides an objective way to identify changes in displacement rates. The displacement rates, their spatial variation, and the timing of accelerations and decelerations can be used to study the physical behavior of a slow-moving slope or for regional hazard assessment by linking the timing of changes in displacement rates to landslide causal and triggering factors.

This book provides a detailed description of how Python can be used to give insight into the flow of groundwater based on analytic solutions. Starting with simple problems to illustrate the basic principles, complexity is added step by step to show how one-dimensional and two-dimensional models of one or two aquifers can be implemented. Steady and transient flow problems are discussed in confined, semi-confined, and unconfined aquifers that may include wells, rivers, and areal recharge. Special consideration is given to coastal aquifers, including the effect of tides and the simulation of interface flow. Application of Python allows for compact and readable code, and quick visualization of the solutions. Python scripts are provided to reproduce all results. The scripts are also available online so that they can be altered to meet site-specific conditions. This book is intended both as training material for the next generation of university students and as a useful resource for practitioners. A primer is included for those who are new to Python or as a refresher for existing users.

Aquifer thermal energy storage (ATES) is an energy efficient technique to provide heating and cooling to buildings by storage of warm and cold water in aquifers. In regions with large demand for ATES, ATES adoption has lead to congestion problems in aquifers. The recovery of thermal energy stored in aquifers can be increased by reducing the distance between wells of the same temperature while safeguarding individual system performance. Although this approach is implemented in practice, the understanding of how this affects both the recovery efficiency and the needed pumping energy is lacking. In this research, the effect of well placement on the performance of individual systems is quantified, and guidelines for planning and design are developed. Results show an increase in thermal recovery efficiency of individual systems when the thermal zones of wells of the same temperature are combined, which is explained by reduced surface area of the thermal zone over which losses occur. The highest increase of the thermal recovery efficiency is found for systems with a small storage volume and long well screens. The relative increase of the thermal recovery efficiency is 12% for average-sized systems with a storage volume of 250,000 m³/year, and 25% for small systems (50,000 m³/year). The optimal distance between wells of the same temperature is 0.5 times the thermal radius, following the trade-off between an increase of the thermal recovery efficiency and the increase in pumping energy. The distance between wells of opposite temperature must be larger than three times the thermal radius to avoid negative interaction.

The estimation of groundwater recharge is of paramount importance to assess the sustainability of groundwater use in aquifers around the world. Estimation of the recharge flux, however, remains notoriously difficult. In this study the application of nonlinear transfer function noise (TFN) models using impulse response functions is explored to simulate groundwater levels and estimate groundwater recharge. A nonlinear root zone model that simulates recharge is developed and implemented in a TFN model and is compared to a more commonly used linear recharge model. An additional novel aspect of this study is the use of an autoregressive-moving-average noise model so that the remaining noise fulfills the statistical conditions to reliably estimate parameter uncertainties and compute the confidence intervals of the recharge estimates. The models are calibrated on groundwater-level data observed at the Wagna hydrological research station in the southeastern part of Austria. The nonlinear model improves the simulation of groundwater levels compared to the linear model. The annual recharge rates estimated with the nonlinear model are comparable to the average seepage rates observed with two lysimeters. The recharges estimates from the nonlinear model are also in reasonably good agreement with the lysimeter data at the smaller timescale of recharge per 10 d. This is an improvement over previous studies that used comparable methods but only reported annual recharge rates. The presented framework requires limited input data (precipitation, potential evaporation, and groundwater levels) and can easily be extended to support applications in different hydrogeological settings than those presented here.

The high landslide risk potential along the steep hillslopes of the Eastern Andes in Ecuador provides challenges for hazard mitigation, especially in areas with hydropower dams and reservoirs. The objective of this study was to characterize, understand, and quantify the mechanisms driving the motions of the Guarumales landslide. This 1.5 km² deep-seated, slow-moving landslide is actively moving and threatening the “Paute Integral” hydroelectric complex. Building on a long time series of measurements of surface displacement, precipitation, and groundwater level fluctuations, we analyzed the role of predisposing conditions and triggering factors on the stability of the landslide. We performed an analysis of the time series of measured groundwater levels and drainage data using transfer functions. The geological interpretation of the landslide was further revised based on twelve new drillings. This demonstrated a locally complex system of colluvium deposits overlying a schist bedrock, reaching up to 100 m. The measured displacement rates were nearly constant at ~50 mm/year over the 18 years of study. However, the measurement accuracy and time resolution were too small to identify possible acceleration or deceleration phases in response to hydro-meteorological forcing. The groundwater and slope drainage data showed a lagged response to rainfall. Finally, we developed a conceptual model of the Guarumales landslide, which we hope will improve our understanding of the many other deep-seated landslides present in the Eastern Andes.

Geogenic arsenic in drinking water is a worldwide problem. For private well owners, testing (e.g., private or government laboratory) is the main method to determine arsenic concentration. However, the temporal variability of arsenic concentrations is not well characterized and it is not clear how often private wells should be tested. To answer this question, three datasets, two new and one publicly available, with temporal arsenic data were utilized: 6370 private wells from New Jersey tested at least twice since 2002, 2174 wells from the USGS NAWQA database, and 391 private wells sampled 14 years apart from Bangladesh. Two arsenic drinking water standards are used for the analysis: 10 µg/L, the WHO guideline and EPA standard or maximum contaminant level (MCL) and 5 µg/L, the New Jersey MCL. A rate of change was determined for each well and these rates were used to predict the temporal change in arsenic for a range of initial arsenic concentrations below an MCL. For each MCL and initial concentration, the probability of exceeding an MCL over time was predicted. Results show that to limit a person to below a 5% chance of drinking water above an MCL, wells that are ½ an MCL and above should be tested every year and wells below ½ an MCL should be tested every 5 years. These results indicate that one test result below an MCL is inadequate to ensure long-term compliance. Future recommendations should account for temporal variability when creating drinking water standards and guidance for private well owners.

Distributed temperature sensing (DTS) systems can be used to estimate the temperature along optic fibers of several kilometers at a sub-meter interval. DTS systems function by shooting laser pulses through a fiber and measuring its backscatter intensity at two distinct wavelengths in the Raman spectrum. The scattering-loss coefficients for these wavelengths are temperature-dependent, so that the temperature along the fiber can be estimated using calibration to fiber sections with a known temperature. A new calibration approach is developed that allows for an estimate of the uncertainty of the estimated temperature, which varies along the fiber and with time. The uncertainty is a result of the noise from the detectors and the uncertainty in the calibrated parameters that relate the backscatter intensity to temperature. Estimation of the confidence interval of the temperature requires an estimate of the distribution of the noise from the detectors and an estimate of the multi-variate distribution of the parameters. Both distributions are propagated with Monte Carlo sampling to approximate the probability density function of the estimated temperature, which is different at each point along the fiber and varies over time. Various summarizing statistics are computed from the approximate probability density function, such as the confidence intervals and the standard uncertainty (the estimated standard deviation) of the estimated temperature. An example is presented to demonstrate the approach and to assess the reasonableness of the estimated confidence intervals. The approach is implemented in the open-source Python package “dtscalibration”.