Control Volume Isogeometric Analysis for groundwater flow modeling in heterogeneous porous media
Introduction
The analysis of flow in heterogeneous porous media is not only important due to water supply or interpretation of hydrogeological properties. Currently, analyses of contaminant transport, remediation or energy exploitation have become very important due to sustainable development and water resources protection. Realistic and reliable transport analysis requires extensive flow simulations including large domains, complex geometries, extensive computational burden with high parallel computing and resolving of many (unfortunately especially fine) heterogeneity scales. Finally, the main input for transport analysis is an accurate velocity field that satisfies the local and global mass balance property and enables continuity of velocity and its derivatives due to the description of a velocity dependent dispersion tensor (Cainelli et al, 2012).
Over the last few decades, many numerical and stochastic methods have been developed to fulfil all the mentioned requirements. Due to scarcity of input data, especially related to the hydraulic conductivity, all flow and transport fields/variables, such as hydraulic conductivity, head, velocity, concentration and/or travel time, have been regarded as Stochastic Random Fields (SRF) (see comprehensive reviews of Dagan, 1989, Gelhar, 1993, Rubin, 2003 and Zhang, 2002). Therefore, many stochastic methods have been developed, such as spectral methods (e..g. Gelhar, 1993), perturbation methods (e.g. Dagan, 1989, Cvetkovic et al., 1992), moment methods (e.g. Hsu et al., 1996), stochastic finite element (e.g. Aguirre and Haghighi, 2003) or probabilistic collocation methods (e.g. Li and Zhang, 2007; Xiu and Hesthaven, 2018) as well as many others. Most of these methods are faced with many computational and implementation problems, especially in the case of high heterogeneity and large anisotropy in irregular domains. Despite the many desired properties of the aforementioned methodologies, the Monte Carlo method is still the most general and reliable stochastic method which at least serves as a verification method of all the mentioned methods. Among others, the Monte Carlo method has been used for flow and transport analysis in weakly and highly heterogeneous porous media by (Bellin et al, 1992, Cvetkovic et al., 1992, 1996, Salandin and Fiorroto, 1998, Zinn and Harvey, 2003, Jankovic et al., 2006, de Dreuzy et al., 2007, Gotovac et al., 2009b, Beaudoin and de Dreuzy, 2013, Ramasomanana et al., 2013, Srzic et al., 2013, Fiori et al., 2015).
The Monte Carlo method first generates as many conditional or unconditional heterogeneity realizations of hydraulic conductivity as possible. Then, the flow equation is solved for all realizations to obtain the velocity field that is the main input in transport simulations of, for example, contaminant transport. The next step is solving the transport equation for all flow realizations and obtaining the plume concentration inside the domain or travel time/flux breakthrough curve at some chosen control planes. Finally, ensemble statistics for all flow and transport variables are analyzed to obtain its statistical moments, correlations or probability density functions – pdf's (Rubin, 2003). However, the computational heart of the Monte Carlo method as a stochastic methodology lies in applying the appropriate numerical techniques for flow and transport simulations in each realization. In this paper, we will focus on numerical flow simulations and derivation of an accurate velocity field needed for reliable transport simulations.
Cainelli et al. (2012) recently analyzed convergence and approximation properties of common numerical methods for flow analysis in heterogeneous porous media. They defined common methods as the widely used finite element (FEM) and finite volume (FVM) methods, which regard hydraulic conductivity as a piecewise element/cell constant field (discontinuous field) representing heterogeneity at some chosen scale.
The classic finite element method (FEM) (e.g. Bellin et al., 1992) divides the domain into the finite elements describing the head as a linear combination of classic Lagrangian basis functions (usually linear or quadratic). The inherent problem of the classic FEM is discontinuity of head gradients along the inter-element boundaries due to the properties of Lagrangian basis functions. This effect, in combination with the common assumption of the discontinuous conductivity field through the Darcy law, produces numerically “unphysical” velocity discontinuities (variations) causing additional artificial numerical dispersion and for instance inaccurate particle tracking in transport simulations. These problems arise with higher heterogeneities due to larger inter-element conductivity contrasts and larger velocity numerical errors. Classic FEM can also violate local mass balance over elements. The first way to overcome classic FEM problems is to employ the mixed finite element (MFE) formulation (e.g. Mose et al. 1994), where velocity is also represented by basis functions involving additional degrees of freedom and computational work in flow simulations. Additionally, some MFE formulations can create ill-posed problems. Another way is to perform postprocessor calculations on FEM velocity fields requiring continuity of velocities between elements and/or mass continuation (Yeh et al., 1981; Cordes and Kinzelbach, 1992).
On the other side, FVM methods, such as MODLFLOW formulation [e.g. Harbaugh, 2005, de Dreuzy et al., 2007; can also be regarded as cell finite difference formulation], satisfy mass balance over the cell boundaries and ensure velocity continuation normal to the cell boundaries. However, MODFLOW generates discontinuities of tangential velocity components between cells. It is possible to prove equivalence between MODFLOW and MFE which uses low order basis functions (Canielli et al., 2012). In addition of classical MODFLOW and FEM, most other existing methods (e.g. Hou and Wu, 1997, Jankovic et al., 2006, Vent and Kitanidis, 1996 or Rasaei and Sahimi, 2009) do not support a true multiresolution approach respecting that hydraulic conductivity can be spanned by many scales.
Gotovac et al. (2009a) presented a such multiresolution approach based on the Adaptive Fup Collocation Method (AFCM) with the following crucial properties: a) correlated log-conductivity (lnK) field is presented as a continuous function in a multiresolution way that can accurately describe all needed scales, b) independently resolves all heterogeneity and head scales, c) yields a continuous velocity field as well as its derivatives needed for application to realistic transport simulations with a velocity dependent dispersion tensor and d) actively controls head and velocity numerical errors reducing numerical dispersion and ensuring convergence, even in highly heterogeneous porous media. The heart of the methodology lies in Fup basis functions (closely related to the B-splines, see Gotovac and Kozulic, 1999 and Supplementary material) in conjunction with the collocation procedure. However, the main drawback was the lack of local and global mass balance due to the properties of the collocation framework, inability to describe the general irregular geometry and computationally expensive head solution to obtain an accurate velocity field without numerical oscillations for high heterogeneity cases. Among others, applications of the AFCM have been shown for the analysis of the flow and transport in heterogeneous porous media relating to the travel time statistics (Gotovac et al., 2009b), concentration statistics (Srzic et al., 2013) and advective transport (Fiori et al., 2015).
Therefore, in this paper, we put desirable AFCM properties to the broader context of isogeometric analysis (IGA). Initially published by Hughes et al. (2005), IGA has attracted enormous attention by designers and modelers because it bridges the gap between Computational Aid Design (CAD) of irregular domains and the numerical description of physical processes on them. The idea is quite simple: use the same type of spline basis functions (usually B-splines or NURBS) to describe both the geometry and solution (head and velocity for our particular groundwater flow problem). IGA has been successfully applied in many fields of solid, structural and fluid mechanics (e.g. Cottrell, Hughes and Bazilevs, 2009). It has been shown that IGA significantly outperforms classical FEM due to higher accuracy and enhanced continuity. Since classic IGA uses the Galerkin variation principle for solution description and transforms the physical domain to the virtual one for geometry description, the apparent similarity between IGA and FEM has resulted in the scientific community widely accepting IGA as new powerful modeling machinery (spline FEM).
In spirit, IGA is closely related to the meshless or mesh-free methodologies due to its use of spline basis functions, which enable many properties not seen in FEM, such as exact geometry description in CAD sense, no classical meshing, more efficient refinement adaptive procedures and multiresolution approach, higher continuity of solution and geometry and usage of higher-order basis functions. In addition to the Galerkin weak formulation, IGA recently presented a collocation approach (e.g. Aurichio et al., 2010; Shillinger et al., 2013; Gomez and De Lorenzis, 2016). Recently, our numerical group developed and applied control/finite volume IGA formulation to the groundwater flow modeling in karst aquifers which are represented by 3-D heterogeneous porous matrix and 1-D laminar/turbulent karst conduits (Malenica et al., 2018).
Efficient numerical modeling using spline basis functions may not always link exclusively with IGA that includes geometry transformations, because everything can be performed only in the physical domain. In that case, spline basis functions which intersect domain boundaries must satisfy the geometry constraints and boundary conditions using for instance the Rvachev solution structure method (see for instance Rvachev et al., 2000, Tsukanov and Shapiro, 2005, Hollig et al., 2003, Apprich et al., 2016, Kozulic and Gotovac, 2017).
Recently, IGA has been used for solving groundwater flow problems such as: filtration through the dam foundation (Shahrbanozadeh et al., 2015), saturated flow (Bekele et al., 2016), unsaturated flow (Nguyen et al., 2014) or flow in fractured media (Hageman and Borst, 2019).
Finally, in this paper, using the framework of Control Volume IsoGeometric Analysis (CV-IGA), the groundwater flow in heterogeneous porous media is then analyzed using two types of spline basis functions, B-splines and Fup basis functions, and compare to the common IGA choices: Galerkin (G-IGA) and Collocation (C-IGA). Except geometry and solution (head and velocity), heterogeneity is also described by Fup/B-spline basis functions introducing CV-IGA as a promising tool for more powerful flow simulations in heterogeneous porous media, which enables the derivation of the velocity field with all the desired mentioned properties. Therefore, this methodology can open new avenues for more robust and realistic simulations of tracer and contaminant transport (Cai et al., 1997, Zheng and Bennett, 2002).
Section snippets
Problem formulation
This section will present brief review of mathematical and physical formulations of flow in heterogeneous porous media. We assume that groundwater flow in heterogeneous porous media can be described by Darcy law:where q is the vector of Darcy velocity, h is the hydraulic head and K is the hydraulic conductivity symmetric tensor of the second order:where Kij is the tensor component which characterizes velocity in i direction under the unit hydraulic gradient in the j
Methodology
This section will provide a detailed description of the methodology employed. There are two basic categories of spline based methods: a) with geometry mapping from the physical space to the parameter space, which is widely coined as IsoGeometric Analysis (IGA), see Hughes et al. (2005) and Cottrell et al. (2007, 2009), and b) without geometry mapping where the solution is directly approximated in the irregular physical domain, see Kozulic and Gotovac, 2017, Tsukanov and Shapiro, 2005 or
Flow problem in 2-D uncorrelated heterogeneous porous media with pumping
This example will show verification of CV-IGA using low-order Fup1 basis functions with respect to the classical methods: MODFLOW and FEM. Verification test is taken from the MODFLOW report (Mehl and Hill, 2005). It is 1260 m long and 1350 m wide, has constant-head boundaries of 10 m and 1.0 m on the left and right sides, respectively. No-flow boundaries extend across the top and bottom (Fig. 4). A pumping well extracts 5.5 m3/s from the center of the domain. The heterogeneity structure was
Conclusions
In this paper, we apply the CV-IGA framework which yields velocity field with all the desired properties that are indispensable for reliable transport simulations:
- –
Consider heterogeneous hydraulic conductivity as a lnK continuous function, resolving all needed spatial heterogeneity scales.
- –
Convergence rate for velocity is optimal (p=n) irrespective to the heterogeneity level if correlated log-conductivity field (16) is defined as a continuous function.
- –
Yields a continuous velocity field (Eq. 30)
Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Acknowledgments
This research is funded by the Croatian Science Foundation (in Croatian Hrvatska zaklada za znanost- HRZZ) through the scientific project “Groundwater flow modeling in karst aquifers”, grant number UIP-2013-11-8103.
This research is partially supported through project KK.01.1.1.02.0027, a project co-financed by the Croatian Government and the European Union through the European Regional Development Fund - the Competitiveness and Cohesion Operational Programme.
References (69)
- et al.
Adaptively refined multi-patch B-splines with enhanced smoothness
Appl. Math. Comput.
(2016) - et al.
On the accuracy of classic numerical schemes for modeling flow in saturated heterogeneous formations
Adv. Water Res.
(2012) - et al.
Fast formation of isogeometric Galerkin matrices by weighted quadrature
Comput. Methods Appl. Math. Eng.
(2017) - et al.
Studies of refinement and continuity in isogeometric structural analysis
Comput. Meth. Appl. Mech. Eng.
(2007) - et al.
PetIGA: a framework for high-performance isogeometric analysis
Comput. Meth. Appl. Mech. Eng.
(2016) - et al.
Isogeometric collocation: Neumann boundary conditions and contact
Comput. Meth. Appl. Mech. Eng.
(2015) - et al.
Algorithms for the implementation of adaptive isogeometric methods using hierarchical B-splines
Appl. Numer. Math.
(2018) - et al.
The variational collocation method
Comput. Meth. Appl. Mech. Eng.
(2016) - et al.
Multi-resolution adaptive modeling of groundwater flow and transport problems
Adv. Water Res.
(2007) - et al.
Adaptive Fup multi-resolution approach to flow and advective transport in highly heterogeneous porous media: Methodology, accuracy and convergence
Adv. Water Res.
(2009)