Control Volume Isogeometric Analysis for groundwater flow modeling in heterogeneous porous media

https://doi.org/10.1016/j.advwatres.2020.103838Get rights and content

Highlights

  • Application of the control volume isogeometric analysis (CV-IGA) to the groundwater flow modeling in heterogeneous porous media.

  • Heterogeneity multiresolution distribution by CV-IGA as a continuous lnK function resolving all resolution scales.

  • Convergence analysis for CV-IGA in comparison to other IGA formulations and classic methods such as MODFLOW and FEM.

  • Obtaining velocity field as a smooth continuous field with optimal convergence rate irrespective to the heterogeneity level.

Abstract

Many important groundwater transport applications require solving the Darcy flow in heterogeneous porous media. Flow simulations, especially in large, highly heterogeneous aquifers, require extensive computational resources, a multiresolution (multiscale) approach to resolve the different heterogeneity scales and an accurate calculation of the velocity field. Common methods, such as finite volumes and elements, assume a discontinuous conductivity field introducing velocity discontinuities along the cell or element interfaces due to using classic discrete operators or Lagrangian basis functions. Over the last decade, the development of isogeometric analysis (IGA) eliminates many of the aforementioned limitations bridging the gap between CAD and numerical analysis. Since classic IGA uses the Galerkin and collocation approach, in this paper, we present a third concept in the form of Control Volume IsoGeometric Analysis (CV-IGA) enabling local and global mass conservation as well as multiresolution description of all heterogeneity scales. Due to the approximation properties of spline basis functions, the velocity field and its derivatives are continuous and are obtained by an optimal convergence rate. The CV-IGA methodology is verified with 2-D numerical and stochastic benchmark flow simulations, including comparisons with classic methods and two other IGA formulations as well as the convergence analysis of the head and velocity fields for different orders of Fup and B-spline basis functions.

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:q=K·hwhere q is the vector of Darcy velocity, h is the hydraulic head and K is the hydraulic conductivity symmetric tensor of the second order:K=[KxxKxyKyxKyy]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)

  • P. Hennig et al.

    Adaptive mesh refinement strategies in isogeometric analysis - a computational comparison

    Comput. Meth. Appl. Mech. Eng.

    (2017)
  • T.Y. Hou et al.

    A multiscale finite element method for elliptic problems in composite materials and porous media

    J. Comput. Phys.

    (1997)
  • T.J. Hughes et al.

    Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement

    Comput. Meth. Appl. Mech. Eng.

    (2005)
  • M.N. Nguyen et al.

    Isogeometric analysis for unsaturated flow problems

    Comput. Geotech.

    (2014)
  • V.P. Nguyen et al.

    Isogeometric analysis: an overview and computer implementation aspects

    Math. Comput. Simul

    (2015)
  • D. Schillinger et al.

    Isogeometric collocation: cost comparison with Galerkin methods and extension to adaptive hierarchical NURBS discretizations

    Comput. Meth. Appl. Mech. Eng.

    (2013)
  • M. Shahrbanozadeh et al.

    Simulation of flow through dam foundation by isogeometric method

    Eng. Sci. Technol. Int. J.

    (2015)
  • A.-V. Vuong et al.

    A hierarchical approach to adaptive local refinement in isogeometric analysis

    Comput. Meth. Appl. Mech. Eng.

    (2011)
  • C.G. Aguirre et al.

    Stochastic finite element analysis of transient unsaturated flow in porous media

    Trans. Am. Soc. Agric. Eng.

    (2003)
  • C. Apprich et al.

    Collocation with WEB–Splines

    Adv. Comput. Math.

    (2016)
  • F. Aurichio et al.

    Isogeometric collocation methods

    Math. Models Methods Appl. Sci.

    (2010)
  • B.R. Baliga et al.

    A control volume finite-element method for two-dimensional fluid flow and heat transfer

    Numer. Heat Transf.

    (1983)
  • Y. Bazargan-Lari

    A pointwise approach for enforcement of essential boundary conditions in the isogeometric analysis

    IJST, Trans. Mech. Eng.

    (2014)
  • Michael Barton et al.

    Optimal quadrature rules for isogeometric analysis

    arXiv: Numer. Anal.

    (2015)
  • A. Beaudoin et al.

    Numerical assessment of 3-D macrodispersion in heterogeneous porous media

    Water Resour. Res.

    (2013)
  • Y.W. Bekele et al.

    Adaptive isogeometric finite element analysis of steady-state groundwater flow

    Int. J. Numer. Anal. Methods Geomech.

    (2016)
  • A. Bellin et al.

    Simulation of dispersion in heterogeneous porous formations: statistics, first-order theories, convergence of computations

    Water Resour. Res.

    (1992)
  • A. Bellin et al.

    Hydro_gen: A spatially distributed random field generator for correlated properties

    Stochastic Hydrol. Hydraulics

    (1996)
  • Z. Cai et al.

    Control volume mixed finite elements

    Comput. Geosci.

    (1997)
  • C. Cordes et al.

    Continuous groundwater velocity fields and path lines in linear, bilinear, and trilinear finite elements

    Water Resour. Res.

    (1992)
  • J.A. Cottrell et al.

    Isogeometric Analysis Toward Integration of CAD and FEA

    (2009)
  • V. Cvetkovic et al.

    A solute flux approach to transport in heterogeneous formations: 2. uncertainty analysis

    Water Resour. Res.

    (1992)
  • V. Cvetkovic et al.

    Analysis of nonlinear effects on tracer migration in heterogeneous aquifers using lagrangian travel time statistics

    Water Resour. Res.

    (1996)
  • D. Dagan

    Flow and Transport in Porous Formations

    (1989)
  • Cited by (0)

    View full text