On proper orthogonal decomposition (POD) based reduced-order modeling of groundwater flow through heterogeneous porous media with point source singularity

Highlights

• Represents local heterogeneity at micro-scale level with an average error variance less than 0.02.

• The FSM shows consistent results at well node locations with changing grid dimensions.

• ROM is computationally efficient and replicates FSM with desired accuracy.

• Both FSM and ROM shows good mesh convergence.

Abstract

Groundwater, being a vital component of the natural water resource system, needs continuous monitoring and dynamic management strategies. That said, we require computationally inexpensive groundwater flow models for repetitive solutions with desirable accuracy under budgetary limitation(s). Natural aquifer systems inherit strong heterogeneity at local scales. In this work, we have proposed ordinary kriging-based sequential algorithm for generating replicates of randomly distributed heterogeneous hydraulic conductivity field (Monte Carlo method-based algorithm) conditioned by field values from sampled locations in an irregular-unstructured grid system. Finite Volume method-based groundwater models often encounter difficulties with the representation of point source/sink terms operating within the domain. In this paper, we have proposed an irregular-unstructured grid Finite Volume discretization technique for overcoming the singularity of point source/sink term to yield a consistent output with different grid dimensions. Furthermore, full-system groundwater models often come with a substantial computational burden. Hence, reduction in model order cuts down the computational expenses (in terms of CPU time and usage) to a significant level. We have also put forth a model order reduction methodology for three different illustrative pumping tests. The proposed framework for the model order reduction projects the governing groundwater flow equation onto a set of identified patterns or orthonormal basis functions, applying the Galerkin Projection method to compute a vector of time-dependent coefficients. We have performed pattern identification by Singular Value Decomposition (SVD) of snapshots of full-system model simulation data at selected time instants within the pumping test time domain. The numerical results of the proposed reduced-order models show a good approximation of the full-system models at a comparatively lesser computational time. The accuracy and efficiency of the models attempt to ensure their potential applicability for identifying groundwater dynamics.

Introduction

Groundwater management framework requires the repetitive use of groundwater flow models for designing a sustainable policy. With time, scientists have developed several mathematical models to estimate the groundwater dynamics of different aquifer systems. Finite Difference (FD) methods-based models such as MODFLOW (McDonald, Harbaugh, 1988, Harbaugh, Banta, Hill, McDonald, 2000, Harbaugh, 2005) have wide applications in groundwater modeling. As FD models work on a simple rectangular grid system, the discretization of governing equations is very straightforward. This makes the FD scheme popular among practitioners. However, the FD scheme fails to describe precisely the boundaries of domains with complex geometry, which significantly affects the accuracy of model prediction near the boundaries. Although grid refinement can be a remedy for the incurred loss in accuracy (Mehl and Hill, 2002), it comes at the cost of increased computational expenses. An improved version of the FD scheme, the Integrated Finite Difference (IFD) method also known as the Control Volume Finite Difference (CVFD) method, has been used as well to analyze fluid flow in porous media (Narasimhan and Witherspoon, 1976). MODFLOW-USG (Panday et al., 2013), a CVFD-based unstructured grid version of MODFLOW, has been developed to solve tightly coupled multiple hydrological flow processes. Other numerical methods like Finite Element (FE) method (Reddy, 2006), Boundary Element (BE) method (Kobayashi and Nishimura, 1992), Finite Volume (FV) method (Leveque, 2002, Mazumder, 2016), and different Mesh-free methods (Liu and Gu, 2005) have proved their relevance in modeling groundwater flow over time. Another attractive numerical method for groundwater modeling is the FE method. The FE-based commercial computer models like FEFLOW (Diersch, 2014) and FEMWATER (Lin et al., 1997) are also available for simulation of saturated/unsaturated subsurface flow and transport. Nevertheless, the accuracy of the FE method suffers due to its non-conservative nature. Nowadays, the FV method is gaining acceptance among emerging numerical methods because it is highly conservative. Moreover, the FV method is also applicable to domains with complex geometries. The FV method integrates the flow-governing equations over the control volume of any shape and conserves the mass, momentum, and the energy entering and exiting a system. FV models for groundwater flow applications have been developed for different regular and irregular-unstructured grid systems (Loudyi, Falconer, Lin, 2006, Dotlić, Vidović, Pokorni, Pušić, Dimkić, 2016).

In groundwater modeling, a well is treated as a singular point source or sink. Moreover, a well encountered in groundwater hydrology, be it of recharging or extraction type, is of vanishingly small radius in the order of 0.1-0.2 metres (Bear, 1979, Pinder, Celia, 2006). The drawdown inside such a small diameter well can be calculated neglecting the water derived from the storage inside the well (Papadopulos and Cooper, 1967). In numerical modeling, a well is specified at the centroid of an element when the governing partial differential equation is discretized on a regular-structured or unstructured grid using FD and CVFD schemes. The problem arises during simulation of real-life situations because the circumference of a well where the pressure gradient is maximum is practically very small compared to the area of the element containing it. Peaceman method defines the well-block pressure in terms of the pressure in the neighbouring elements. It is depicted to be numerically equal to the actual pressure at an effective radius of 0.2Δx for regular square grid (Peaceman, 1977) and $0.14{({\Delta x}^2+{\Delta y}^2)}^{\frac{1}{2}}$ for regular rectangular grid (Peaceman, 1983) systems. However, local grid refinement around a well can be considered to bring down the effective radius very close to the actual value. On the other hand, finer spatial discretization imposes several restrictions on temporal discretization to maintain stability of the numerical scheme. This eventually increases the computational expenses turning it into an impractical attempt. Furthermore, when irregular-unstructured grid systems are considered, the above considerations do not function as grid dimensions vary for each element, thereby separate well models become a necessity. Different well models for standard FEM, control volume FEM, and mixed FEM have been developed by Chen and Zhang (2009). The well model derived for control volume FEM has been applied for CVFD scheme in MODFLOW-USG (Panday et al., 2013) which is approximately equivalent to Peaceman method applied for structured grids in previous versions of MODFLOW. Regardless, bringing down the effective well radius close to the actual radius value still remains an unanswered question.

The application of any groundwater flow model requires the identification of aquifer parameters such as the hydraulic conductivity field and the storage coefficients. The primary means of aquifer characterization are on-field pumping tests and numerical simulations. As pumping tests are complex and expensive, alternative economical methods based on particle-size distribution (Rogiers et al., 2012) and Electrical Resistivity Tomography (ERT) data (Yu and Wu, 2006) have been used to estimate field hydraulic conductivity. One of the significant challenges faced in modeling groundwater flow through porous media is the mathematical replication of parameter heterogeneity inherent to transport media. Previously, to avoid complexities, the concept of equivalent porous medium was prevalent where the effects of heterogeneity were lumped into a single effective parameter which represented the medium. Yeh (1986) has modeled the field heterogeneity by dividing the flow regime into several sub-regions where a constant hydraulic conductivity characterizes each region. Kriging is an important tool used in groundwater hydrology for distributing field information throughout the domain from limited data points. The kriging estimates are highly influenced by the uncertainty of the semivariogram as it is determined on the basis of available sampled data. The uncertainty of the model semivariogram can be quantified using random fields which is very essential for stochastic modeling (Lin and Chen, 2005). Xu et al. (2017) has shown the successful application of kriging-based interpolation technique to determine spatial distribution of hydraulic conductivity obtained from resistivity and grain-size data at sampled locations.

In reality, aquifer parameters not only vary among different layers, the variation within a particular layer is also of much significance. The heterogeneity of a hydraulic conductivity field is stochastically best described by a log-normal probability density function (Freeze, 1975). A stochastic model of porous media considers hydraulic conductivity as a random variable (Greenkorn and Kessler, 1969), and the intrinsic heterogeneity as a Space Random Function (SRF) correlated in space and time (Bellin and Rubin, 1996). Monte Carlo method-based sequential Gaussian simulation is considered an efficient algorithm for the generation of random fields. Application of the Monte Carlo method to such SRF field generation simplifies the intricacies associated with the replication of spatial variability of field heterogeneity. However, a significant limitation of the Monte Carlo method is that it converges slowly as $\mathcal{O}\left[\frac{1}{\sqrt{N}}\right]$ where N is the number of independent Monte Carlo realizations. In order to obtain nearly accurate results, we require a large number of Monte Carlo replicates of the random field, which eventually increases the computation time of the model. Many effective commercial softwares are available for random field generation based on sequential Gaussian simulation. Random Surface Generation (RSG) is a Moving Average method-based software used for generating Gaussian random surfaces to account for the topographic height variations over a surface (Bergström et al., 2007). Similarly, Random Finite Element Method (RFEM) is a FEM-based geotechnical modeling software which generates Gaussian random fields for introducing random behavior of geologic formations (Smith et al., 2013).

On the other hand, exhibiting the characteristics of a hydraulic conductivity field requires the generated field realization to be conditioned by field measurements. Conditional simulation of heterogeneous porous formations based on measured hydraulic conductivity or transmissivity data at few points of the formation have been successfully attempted (Delhomme, 1979, Dagan, 1982, Clifton, Neuman, 1982, Tsai, 2006). TBSIM, a computer program for conditional simulation of Gaussian random fields (Emery and Lantuéjoul, 2006) is also available for simulating random fields conditioned by available data at selected locations. By far, Stanford Geostatistical Modeling Software (SGeMS) (Remy et al., 2009) is the most productive conditional Gaussian random field generator used to generate random fields replicating the properties of heterogeneous porous media. Although the aquifer is intrinsically heterogeneous, in continuum approach to flow through porous media, the parameters are considered to be homogeneous within a small volume known as Representative Elementary Volume (REV) (Bear, 1979, Pinder, Celia, 2006). This implies that natural heterogeneous porous media is homogeneous at microscopic scale. Hence, accurate mathematical formulations should consider a suitable REV for the conditioning data while generating fields that characterize natural geologic formations.

The accuracy of a groundwater model highly depends on the discretization of the spatial and the temporal domain. Finer the discretization, more accurate are the model predictions. This increases the computational expenses making it an inefficient one in the case of high-dimensional models. Therefore, the focus should be on the development of a model that optimizes both the efficiency parameters - accuracy and computation time. Scientists of various genre have acclaimed the model order reduction methodology based on data-driven pattern identification techniques in different fields of scientific study over the past few decades (Sirovich, 1987, Park, Cho, 1996, Park, Chung, Lee, 1999, Haasdonk, Ohlberger, 2011, Hasenauer, Löhning, Khammash, Allgöwer, 2012, Dehghan, Abbaszadeh, 2018, Abbaszadeh, Dehghan, 2020). The method projects the governing equation onto a reduced sub-space of identified patterns, thereby cutting down the number of algebraic equations to be solved for each temporal iteration. Vermeulen et al. (2004) made the first attempt towards low-dimensional modeling of groundwater flow systems. Since then, model reduction technique has been applied to modeling of groundwater flow through both confined (Pasetto, Guadagnini, Putti, 2011, Pasetto, Putti, Yeh, 2013, Boyce, Yeh, 2014) and unconfined aquifers (Boyce, Nishikawa, Yeh, 2015, Stanko, Boyce, Yeh, 2016), groundwater contamination problems (Ilati, Dehghan, 2016, Kani, Elsheikh, 2019) and groundwater optimization and management studies (McPhee, Yeh, 2008, Baú, 2012). In reduced-order modeling, snapshot selection is the key to pattern recognition. Hence, the temporal distribution of snapshots for groundwater flow systems (Siade et al., 2010) should be chosen wisely. The applications of model order reduction for characterizing aquifer dynamics are limited to FD and FE methods till date. As the FV method is a highly useful tool for PDE discretization, reducing the degrees of freedom involved in the classical approach will help in minimizing the computational complexities.

In this paper, we have presented an FV-based computationally inexpensive mathematical model for evaluating the dynamics of groundwater flow through a heterogeneous confined aquifer subjected to singular point source/sink. In the first step, we have developed a fast, effective and accurate algorithm for generating spatially-correlated and randomly distributed heterogeneous hydraulic conductivity field based on available field values from sampled locations by integrating the concept of REV into standard sequential Gaussian simulation technique. We have developed the field for an irregular-unstructured grid system to facilitate the application of advanced numerical tools. Secondly, we have presented an FV formulation to estimate the groundwater response through the generated random heterogeneous field subjected to singularity of point source/sink. In the proposed methodology, we have considered the actual radius of the well which is independent of the grid dimensions. We have carried out a separate study of the proposed FV-based full-system model through a homogeneous confined aquifer to compare the performance of the model with MODFLOW simulation results. Finally, we have performed POD-based reduced order modeling to minimize the computational expenses in terms of CPU time and usage acquired by the proposed FV-based full-system model.