Subsurface Flow Model Research Articles

An efficient approach of meshless node placement in three-dimensional subsurface flow modeling

Meshless methods are valuable tools for real-time subsurface flow modeling and are particularly beneficial for adaptively adjusting node positions when incorporating supplementary measurements. However, an algorithm is required for automatically generating high-quality meshless nodes to adapt to irregularly shaped aquifers and arbitrarily distributed wells. This paper introduces a three-dimensional adaptive meshless node placement technique based on advancing front nodes and proposes a sigmoid exclusion circle to ensure adaptability to specific positions of interest. It is crucial for groundwater simulations to have specified computational points at wellbore locations. The generated meshless nodes were used as computational points in the generalized finite difference method. The node quality was assessed by comparing the numerical solutions to the analytical solution, achieving errors of the order of 10–12. The application scenarios of irregularly shaped aquifers demonstrated the high efficiency of the three-dimensional node placement in generating up to 10,000 nodes in <5 s. The effectiveness of the algorithm was verified using three-dimensional irregular computational domains. The algorithm could be easily adapted to two-dimensional problems, and an irregular boundary example is presented. The validation and application cases highlighted the versatility of the proposed method and established its potential for real-world hydrogeological applications.

Learning CO[formula omitted] plume migration in faulted reservoirs with Graph Neural Networks

Deep-learning-based surrogate models provide an efficient complement to numerical simulations for subsurface flow problems such as CO2 geological storage. Accurately capturing the impact of faults on CO2 plume migration remains a challenge for many existing deep learning surrogate models based on Convolutional Neural Networks (CNNs) or Neural Operators. We address this challenge with a graph-based neural model leveraging recent developments in the field of Graph Neural Networks (GNNs). Our model combines graph-based convolution Long-Short-Term-Memory (GConvLSTM) with a one-step GNN model, MeshGraphNet (MGN), to operate on complex unstructured meshes and limit temporal error accumulation. We demonstrate that our approach can accurately predict the temporal evolution of gas saturation and pore pressure in a synthetic reservoir with impermeable faults. Our results exhibit a better accuracy and a reduced temporal error accumulation compared to the standard MGN model. We also show the excellent generalizability of our algorithm to mesh configurations, boundary conditions, and heterogeneous permeability fields not included in the training set. This work highlights the potential of GNN-based methods to accurately and rapidly model subsurface flow with complex faults and fractures.

Development and Application of a New Open-Source Integrated Surface–Subsurface Flow Model in Plain Farmland

Accurately characterizing rainfall runoff processes in plain farmland, especially at the plot scale with significant micro-topographic features, has presented challenges. Integrated surface–subsurface flow models with high-precision surface flow modules are appropriate tools, yet open-source versions are rare. To address this gap, we proposed an open-source integrated surface–subsurface flow model called the FullSWOF-Plain model, in which the one-dimensional subsurface module Hydrus-1D was integrated with a modified two-dimensional surface water flow module (FullSWOF-2D) using the sequential head method. Various experimental scenarios were simulated to validate the model’s performance, including two outflow cases (i.e., 1D and 2D) without infiltration, a classical one-dimensional infiltration case, and two typical rainfall events at the experimental field. The results demonstrate the accuracy of this proposed model, with the Nash–Sutcliffe efficiency (NSE) of the simulated discharge exceeding 0.90 in the experimental field case and the root mean squared error (RMSE) values for soil moisture at five depths consistently below 0.03 cm3/cm3. However, we observed a lag in the simulated response time of soil moisture due to the neglect of preferential flow. The micro-topography significantly influenced ponding time and ponding areas. Lower local terrain normally experienced earlier surface ponding. Scattered surface ponding water first occurred in the ditch and followed in the relatively low areas in the main field. The concentration process of surface runoff exhibited hierarchical characteristics, with the drainage ditch contributing the most discharge initially, followed by the connection of scattered puddles in the main field, draining excess surface water to the ditch through rills. This quantitative study sheds light on the impact of micro-topography on surface runoff in plain farmland areas.

Subsurface flow model for estimating 2D wetting pattern in drip irrigation

AbstractNumerical model using 2D Richard’s equation in cylindrical polar coordinate system was developed to assess the soil wetting patterns in a drip irrigation system. A fully implicit finite‐difference scheme was used to solve equations with suitable initial and boundary conditions. The developed model was validated with an experimental result available in the literature, which shows a high level of concordance with low values of RMSE and R2 greater than 0.93 for all discharge rates. Furthermore, the versatility and applicability of the model were demonstrated for complex subsurface conditions considering homogeneous and heterogeneous soil profiles. Water ponding on the surface was observed in clay loam soil, contrasting with deep percolation in sandy loam soil. This underscores the significance of adjusting both discharge rate and irrigation frequency in accordance with prevailing soil conditions.

Surface-subsurface filtration transport with seawater intrusion: multidomain mixed variational evolution problems

Surface-subsurface filtration transport with seawater intrusion phenomena are formulated and variationally analyzed, as coupled multimedia mixed pairs of free boundary interface problems. Physically, multidomain subsurface mixed velocity-pressure fractional Darcian flow models coupled with surface evolution Stokesian mixed flows are considered. Specifically, two-phase air-fresh water above the sea level and fresh water-seawater characterizations are considered. Internal boundary synchronizing transmission conditions of multidomain nonoverlapping decompositions are modeled in terms of variational Lagrangian dual subpotential maximal monotone inclusions. Similarly, filtration transport coupling interface transmission constraints are implemented by mass flux-velocity-pressure Lagrange dual multipliers as solutions of subpotential subdifferential equations.

Identification of aquifer heterogeneity through inverse methods

The paper underlines the contributions of Ghislain de Marsily (GdM) to the identification of aquifers heterogeneity using inverse methods mainly for modeling subsurface flow. Inverse methods require an objective function to express the goodness of fit of the chosen model, a parameterization to describe the spatial distribution of model parameters, and a minimization algorithm. The resulting inverse problem, which consists in seeking model parameters’ values that render model outputs close to the observations, is usually unstable. GdM developed seminal ideas for the two key inversion issues that are: to stabilize the inverse problem through regularization, and to parameterize it to reproduce the natural heterogeneity of the subsurface with a limited number of parameters. GdM conducted pioneering works that are the basis of current parameterization methods relying upon adaptive zonation and/or interpolation based on pilot points. We take here the opportunity to highlight the GdM’s contributions inspiring currently used techniques.

Gaussian process regression and conditional Karhunen-Loève models for data assimilation in inverse problems

We present a model inversion algorithm, CKLEMAP, for data assimilation and parameter estimation in partial differential equation models of physical systems with spatially heterogeneous parameter fields. These fields are approximated using low-dimensional conditional Karhunen-Loève expansions (CKLEs), which are constructed using Gaussian process regression (GPR) models of these fields trained on the parameters' measurements. We then assimilate measurements of the state of the system and compute the maximum a posteriori (MAP) estimate of the CKLE coefficients by solving a nonlinear least-squares problem. When solving this optimization problem, we efficiently compute the Jacobian of the vector objective by exploiting the sparsity structure of the linear system of equations associated with the forward solution of the physics problem.The CKLEMAP method provides better scalability compared to the standard MAP method. In the MAP method, the number of unknowns to be estimated is equal to the number of elements in the numerical forward model. On the other hand, in CKLEMAP, the number of unknowns (CKLE coefficients) is controlled by the smoothness of the parameter field and the number of measurements, and is generally much smaller than the number of discretization nodes, which leads to a significant reduction of computational cost with respect to the standard MAP method. To show this advantage in scalability, we apply CKLEMAP to estimate the transmissivity field in a two-dimensional steady-state subsurface flow model of the Hanford Site by assimilating synthetic measurements of transmissivity and hydraulic head. We find that the execution time of CKLEMAP scales nearly linearly as N1.33, where N is the number of discretization nodes, while the execution time of standard MAP scales as N2.91. The CKLEMAP method improved execution time without sacrificing accuracy when compared to the standard MAP method.

Development of inter-grid-cell lateral unsaturated and saturated flow model in the E3SM Land Model (v2.0)

Abstract. The lateral transport of water in the subsurface is important in modulating terrestrial water energy distribution. Although a few land surface models have recently included lateral saturated flow within and across grid cells, it is not a default configuration in the Climate Model Intercomparison Project version 6 experiments. In this work, we developed the lateral subsurface flow model within both unsaturated and saturated zones in the Energy Exascale Earth System Model (E3SM) Land Model version 2 (ELMv2.0). The new model, called ELMlat, was benchmarked against PFLOTRAN, a 3D subsurface flow and transport model, for three idealized hillslopes that included a convergent hillslope, divergent hillslope, and tilted V-shaped hillslope with variably saturated initial conditions. ELMlat showed comparable performance against PFLOTRAN in terms of capturing the dynamics of soil moisture and groundwater table for the three benchmark hillslope problems. Specifically, the mean absolute errors (MAEs) of the soil moisture in the top 10 layers between ELMlat and PFLOTRAN were within 1 %±3 %, and the MAEs of water table depth were within ±0.2 m. Next, ELMlat was applied to the Little Washita experimental watershed to assess its prediction of groundwater table, soil moisture, and soil temperature. The spatial pattern of simulated groundwater table depth agreed well with the global groundwater table benchmark dataset generated from a global model calibrated with long-term observations. The effects of lateral groundwater flow on the energy flux partitioning were more prominent in lowland areas with shallower groundwater tables, where the difference in simulated annual surface soil temperature could reach 0.3–0.4 ∘C between ELMv2.0 and ELMlat. Incorporating lateral subsurface flow in ELM improves the representation of the subsurface hydrology, which will provide a good basis for future large-scale applications.

A practical approach for numerical modeling of a complex and data-limited hydrological system

Abstract By their nature, wetlands represent an ecosystem base for many concurrent heterogeneous interactions, where the mission of numerical modeling requires a wide range of consistent and reliable datasets from various sources, spatially and temporally. Such a mission usually collides with the existence of tremendous missing in time-series dataset(s). In this context, MIKE SHE was used to construct an integrated surface–subsurface flow model for the Paya Indah wetland in Malaysia where huge gaps exist in the historical datasets of water level and flow rate. To calibrate and validate the model to a satisfactory level, a tri-criteria simulation approach was applied to overcome the occasional missing values in these datasets. This goal was accomplished by calibrating the surface water level and channel flow while simultaneously matching the steady-state subsurface portion of the system wherever water table depth data allowed. Quantitatively, the integrated model scored the highest values of R (0.765 − 0.927) and CE (0.748 − 0.828) during the validation. However, large RMSE values were calculated for the flow rate during calibration at SWL2 (outlet; 0.766) and during validation at Langat River (0.780). This bias was attributed to the low or occasional absence of variation in the historical time-series datasets necessary for the simulation process.

Data assimilation for real-time subsurface flow modeling with dynamically adaptive meshless node adjustments

Data assimilation for real-time subsurface flow modeling with dynamically adaptive meshless node adjustments

A preliminary development of a coupled surface and subsurface flow model for swale system

This study presents the preliminary stage of the development of a conjunctive surface-subsurface model for simulating the flow in a swale. The surface flow is modelled by one-dimensional dynamic wave equation. To improve accuracy, a third-order numerical scheme is used to solve the advection terms in the dynamic wave model. Meanwhile, the subsurface flow is modelled as one-dimensional vertical flow through nondeformable porous media without air compression effect. At preliminary stage, the surface and subsurface model are verified separately. The surface model is verified against the experimental data of surface runoff from a simulated rainfall. The surface model showed promising performance in terms of the reproduction of surface hydrograph. The overall discrepancy between the experimental and numerical model result in the reproduction of the hydrograph discharge is around 14%. For the subsurface flow model, the vertical moisture profile of the soil is verified against Philip’s analytical solution. The moisture profile obtained from the numerical model shows 100% agreement with Philip’s solution.

On Groundwater Recharge in Variably Saturated Subsurface Flow Models

AbstractGroundwater models that simulate only saturated flow use groundwater recharge as an input parameter. In contrast, variably saturated subsurface flow models, including integrated surface and subsurface hydrologic models, can jointly simulate the movement of water in the saturated and unsaturated zones. Instead of recharge, they require climate data such as precipitation and potential evapotranspiration. Given that the latter models represent hydrological processes operating throughout the unsaturated zone and at the water table, one might expect that recharge can be readily extracted from them. In this paper, we demonstrate that it is not the case. When the commonly used definitions of groundwater recharge are implemented in variably saturated subsurface flow models, they do not yield meaningful results. Above all, the problems occur because of the storage dynamics in the capillary fringe above the water table. Despite this difficulty, variably saturated subsurface flow models can provide the information required for water resources management directly.

Deep ensemble learning for high-dimensional subsurface fluid flow modeling

The accuracy of Deep Learning (DL) algorithms can be improved by combining several deep learners into an ensemble. This avoids the continuous endeavor required to adjust the architecture of individual networks or the nature of the propagation. This study investigates prediction improvements possible using Deep Ensemble Learning (DEL) to determine four distinct multiscale basis functions in the mixed Generalized Multiscale Finite Element Method (GMsFEM), involving the permeability field as the only input. 376,250 samples were initially generated, filtered down to 367,811 after data pre-processing. A standard Convolutional Neural Network (CNN) named SkiplessCNN and three skip connection-based CNNs named FirstSkipCNN, MidSkipCNN, and DualSkipCNN were developed for the base learners. For each basis function, these four CNNs were combined into an ensemble model using linear regression and ridge regression, separately, as part of the stacking technique. A comparison of the coefficient of determination (R2) and Mean Squared Error (MSE) confirms the effectiveness of all three skip connections in enhancing the performance of the standard CNN, with DualSkip being the most effective among them. Additionally, as evaluated on the testing subset, the combined models meaningfully outperform the individual models for all basis functions. The case that applies linear regression delivers R2 ranging from 0.8456 to 0.9191 and MSE ranging from 0.0092 to 0.0369. The ridge regression case achieves marginally better predictions with R2 ranging from 0.8539 to 0.922, and MSE ranging from 0.009 to 0.0349 because its solution involves more evenly distributed weights.

Quantifying the impact of the urban karst on infiltrated stormwater

Abstract Urbanization alters the flow regime of streams, including increasing the frequency and magnitude of storm flows, along with reducing baseflows. An increasingly common management strategy is stormwater infiltration, which is thought to reduce surface runoff and recharge groundwater and thus restore lost baseflows to streams. Recent research has pointed to considerable uncertainty on the fate of infiltrated stormwater, particularly due to the presence of human-made underground infrastructure – e.g. sewer and water supply pipes and telecommunication cables. Such infrastructure is commonly housed in trenches partially filled with highly permeable material which can cause urban karst like flow conditions. We used a dynamic subsurface flow model (HYDRUS-3D) to predict the impact of the urban karst on the fate of infiltrated stormwater. The model was constructed with the presence of a sewer pipe situated between an infiltration basin and a stream. The model predicted that the impact of the urban karst on infiltrated stormwater increases with higher groundwater levels, and greater contrast between hydraulic conductivity of regional soil and gravel which surrounds the sewer pipe. Results suggest that it is important to consider the impact of the urban karst in cases where the goal of stormwater infiltration is baseflow restoration.

The relationship between fingering flow fraction and water flux in unsaturated soil at the laboratory scale

Preferential flow is a widespread process in unsaturated soil and it is important to incorporate preferential flow in mathematical and numerical models to accurately predict subsurface flow and solute transport. In this study we conduct experiments to relate the fraction of fingering flow (f) with water flux (qs) and average active water saturation (Se∗) which is a key constitutive relation in modeling subsurface preferential flow. Infiltration experiments were conducted in a two-dimensional transparent slab chamber, with different soil grain sizes and initial water saturation. Available experimental data from the literature are also compiled and compared with our new data. Results show that the f-qs relationship is independent on grain size of the soil and whether the sand is initially dry or wet. The exponent of the f-qs relationship (a0) varies in a small range (0.39–0.48) and the exponent of the f-Se∗ relationship (γ) can be predicted from a0 and soil grain size. The Richards equation was combined with the f-qs relationship to model the preferential flow experiments. The results show that the constitutive relationship successfully predicts preferential flow velocity in unsaturated soil.

Impact of parameter updates on soil moisture assimilation in a 3D heterogeneous hillslope model

Abstract. Variably saturated subsurface flow models require knowledge of the soil hydraulic parameters. However, the determination of these parameters in heterogeneous soils is not easily feasible and subject to large uncertainties. As the modeled soil moisture is very sensitive to these parameters, especially the saturated hydraulic conductivity, porosity, and the parameters describing the retention and relative permeability functions, it is likewise highly uncertain. Data assimilation can be used to handle and reduce both the state and parameter uncertainty. In this work, we apply the ensemble Kalman filter (EnKF) to a three-dimensional heterogeneous hillslope model and investigate the influence of updating the different soil hydraulic parameters on the accuracy of the estimated soil moisture. We further examine the usage of a simplified layered soil structure instead of the fully resolved heterogeneous soil structure in the ensemble. It is shown that the best estimates are obtained when performing a joint update of porosity and the van Genuchten parameters and (optionally) the saturated hydraulic conductivity. The usage of a simplified soil structure gave decent estimates of spatially averaged soil moisture in combination with parameter updates but led to a failure of the EnKF and very poor soil moisture estimates at non-observed locations.

Mesh-independent multidimensional coupling of surface and subsurface water flow models

Abstract We present a new approach for multidimensional surface–subsurface flow coupling. The method does not require mesh-to-canal refinement or alignment and is based on numerical characteristics of the intersection of one-dimensional hydraulic object beds and the surface faces of the 3D tetrahedral mesh. The method is validated using widely recognized tilted v-catchment with subsurface and Borden catchment benchmarks.

Upscaling Porous Media Using Neural Networks: A Deep Learning Approach to Homogenization and Averaging

In recent years machine learning algorithms have been gaining momentum in resolving subsurface flow issues related to hydrocarbon flows, Carbon capture utilization and storage, hydrogen storage, geothermal flows, and enhanced oil recovery. This paper presents and attempts to solve subsurface flow problem using neural upscaling method. The neural upscaling method, described in the present work, is a machine learning approach to calculate effective properties in each grid block for subsurface flow modeling. This method is intended to be more accurate than traditional analytical upscaling methods (which are only accurate for layered or homogeneous media) and numerical upscaling methods (which are more accurate for heterogeneous media but involve higher computational cost and are dependent on boundary conditions). The neural upscaling method is based on learning from a large number of geological realizations, which allows it to account for uncertainty in geology. It is also computationally fast and accurate. The method is demonstrated through a series of 2D test cases, and its accuracy is compared to that of analytical and numerical upscaling methods.

Keynote lecture. The debris flow event of 29 October 2018 in the Rio Rotiano (Italy) and its challenges for the mathematical and numerical modelling

The debris flow that interested the Rio Rotiano, a creek located in the Province of Trento (Italy) on 29 October 2018, is an event that, because of its peculiar features, presents formidable challenges in terms of mathematical and numerical modelling. Here we present some results of our research, which required both physical and mathematical modelling, aimed at developing a numerical tool necessary to face the back analysis of the event and the validation of the planned protection works. Exploiting some results of the laboratory tests and coupling an advanced description of the debris flow dynamics (TRENT2D model) over fixed and mobile bed with a sub-surface flow model, we obtained a model that can be defined as a mobile-bed rainfall-runoff model, where the runoff is composed by both water and sediments and represents a new paradigm in the field of debris flow simulations. First applications give promising results but some further developments are required before completing the back analysis of the Rio Rotiano event.

Impact of beach face slope variation on saltwater intrusion dynamics in unconfined aquifer under tidal boundary condition

Density-dependent saltwater flow in coastal aquifers varies with the seaside tidal and beach slope conditions. This paper presents a set of laboratory-scale saltwater intrusion experiments under different beach slopes (15°, 20°, 25°, and 30°) for unconfined aquifer conditions. Experimental porewater pressure measurements were utilized for quantitative analysis. A new G channel-based image analysis technique is used for experimental image analysis and freshwater–saltwater interface identification. Experimental results were numerically validated using the two-dimensional Finite Element subsurface FLOW (FEFLOW) model. The time-varying analysis revealed that saltwater intrusion occurs rapidly on flatter beach slopes. In the case of a steeper beach slope, saltwater intruded less. Steeper slopes take more time to reach a quasi-steady condition. Tidal oscillations alter hydraulic gradient changes across the beach slope. This gradient change generated clockwise circulating saltwater flow within porous media in the inter-tidal zone. This circulating flow resulted in the formation of the upper saline plume (USP). The USP expanded with time and moved in a downward direction. Finally, a deformed elliptic-shaped USP was observed under quasi-steady-state conditions. Submarine groundwater discharge (SGD) pathways move the saltwater region through the intermediate zone of a saltwater wedge and USP on varying beach slopes. It is evident that SGD particles move along the saltwater–freshwater interface (SWI) zone and rise upward (up to the intersection point). The tracer experiments were started after attaining the quasi-steady state condition and continued till the tracer reached the intersection point of saltwater level and sloping beach face. The experimental data sets can be used as benchmark test cases.