## 1. Introduction

Fluid flow in porous media has been extensively studied, both theoretically and computationally, because of its broad applications in different branches of science and engineering. The most popular model of flow of an incompressible fluid in rigid porous media is the Darcy model, which is based on the assumption that the domain contains only one pore-network. Due to the restricting assumptions in the classical Darcy model [Rajagopal, 2007; Nakshatrala and Rajagopal, 2011; Chang et al., 2017], its application has been limited and several modifications and alternative models have been proposed that predict more realistic flow behaviors. In particular, due to the complexity of the pore-structure in many geo-materials such as shale, many studies have focused on developing mathematical models and computational frameworks that consider the presence of two (or more) dominant pore-networks exhibiting different hydro-mechanical properties. Some of the recent studies on multiple pore-networks include [Borja and Koliji, 2009; Choo et al., 2015].

The mathematical models pertaining to the flow in porous media with multiple pore-networks are complex and involve numerous field variables. It is not always possible to derive analytical solutions to these mathematical models, and one has to resort to numerical solutions for realistic problems. Different approaches are available for developing formulations for multi-field mathematical models. Mixed finite element formulations, which offer the flexibility of using different approximations for different field variables, are particularly attractive for multi-field problems. Accurate numerical solutions have been obtained using mixed finite element for various porous media models; for example, see [Masud and Hughes, 2002; Badia and Codina, 2010; Nakshatrala et al., 2006; Nakshatrala and Rajagopal, 2011; Choo and Borja, 2015]. Moreover, many of the mathematical models pertaining to the multiple pore-networks, and in particular, the mathematical model considered in this paper, cannot be written in terms of a single-field variable. Although mixed methods are considered a powerful tool, especially for modeling flow problems in porous media, they suffer from some restrictions. To obtain stable and convergent solutions, a mixed formulation should satisfy the Ladyzhenskaya-Babuška-Brezzi (LBB) stability condition [Babuška, 1973; Brezzi and Fortin, 1991]

. Numerical instability of the solution and probable spurious oscillations in the profile of unknown variables are the main consequences of the violation of this condition. Such drawbacks are observed in many of the existing formulations and highlight the need for developing more robust computational frameworks. In order to resolve numerical instabilities resulting from violation of the LBB condition, computational approaches are divided broadly into two classes

[Franca and Hughes, 1988]: those that satisfy the LBB condition and those that circumvent it.In the former approach, elements are developed by placing restrictions on the interpolation spaces so as to satisfy the LBB condition under the classical mixed (Galerkin) formulation. Such elements are collectively referred to as the H(div) elements [Brezzi and Fortin, 1991; Brezzi et al., 2008]. Two popular works of this type are Raviart-Thomas (RT) spaces [Raviart and Thomas, 1977], and Brezzi-Douglas-Marini (BDM) spaces [Brezzi et al., 1985, 1987]. The class of stabilized methods, which falls under the latter approach, is an attractive way of circumventing the LBB condition. In a stabilized formulation, stabilization terms are augmented to the classical mixed formulation to avoid a saddle-point problem as well as mathematical instabilities [Hughes et al., 2004]. Various stabilized formulations have been published for flow problems (e.g., see [Badia and Codina, 2009; Brooks and Hughes, 1982; Hughes et al., 2000; Turner et al., 2009]) and for flow problems in porous media, in particular, (e.g., see [Masud and Hughes, 2002; Badia and Codina, 2010; Nakshatrala et al., 2006; Choo and Borja, 2015]).

Herein, we develop a stabilized mixed formulation of the double porosity/permeability model proposed by [Nakshatrala et al., 2016]. The stabilization terms and the stabilization parameter are derived in a mathematically consistent manner by appealing to the variational multiscale formalism [Hughes, 1995]. It is noteworthy that the nodal-based equal-order interpolation for all the field variables is stable under the proposed stabilized mixed formulation. Such a feature for interpolations is particularly desirable for studies in porous media for two reasons. The obvious reason is that the equal-order interpolation is computationally the most convenient. The second reason is that, in many porous media applications, the flow and transport equations are coupled (Section 9 of this paper deals with such a coupled problem). But many existing formulations (including the stabilized formulations) produce non-physical negative solutions for the transport equations (i.e., a negative value for concentration fields), especially when the diffusion/dispersion is anisotropic [Nagarajan and Nakshatrala, 2011]. The known robust non-negative finite element based formulations for the transport equations are nodal-based (e.g., refer to [Nagarajan and Nakshatrala, 2011; Mudunuru and Nakshatrala, 2016]). By choosing nodal-based unknowns even for the flow problem, one can avoid projections from nodal to non-nodal interpolation spaces and vice-versa.

To determine whether a computational framework is robust, systematic convergence and error
analyses are required.
To this end, we first perform a mathematically
rigorous stability analysis of the proposed
stabilized mixed formulation. Since the proposed
formulation is residual-based, consistency is
shown quite easily. We also present patch tests
and representative numerical results to
show that the obtained numerical results are stable.
After establishing the stability of the proposed
formulation, we perform a thorough accuracy assessment
of the approximations by estimating the error
associated with the numerical solutions.
Specifically, we perform both *a priori*
and a posteriori error estimations,
which individually serve different purposes
[Babuška et al., 2010].
A posteriori error estimations monitor
different forms of the error in the numerical
solution [Becker and Rannacher, 2001; Babuška and Strouboulis, 2001]
and using the computed approximate solution, they provide
an estimate of the form ,
where is the solution, is the finite element
solution for a mesh with mesh size ,
denotes an appropriate norm, and is a
constant (real) number.
On the other hand, a priori error
estimations provide us with the order of
convergence of a given finite element
method [Ainsworth and Oden, 1997].

[Shabouei and Nakshatrala, 2016] have shown that
porous media models such as those defined by the Darcy
and Darcy-Brinkman equations satisfy certain mechanics-based
properties, and they have utilized these properties to construct
solution verification procedures.
Recently, [Nakshatrala et al., 2016]
have shown that the double porosity/permeability model
also enjoys properties with strong mechanics underpinning.
These include
the minimum dissipation theorem and a reciprocal
relation. Herein, we utilize these mechanics-based
properties to construct *a posteriori* solution
verification procedures to assess the accuracy of
numerical solutions obtained under the proposed
formulation for the double porosity/permeability
model.

Another type of numerical instability, known as Gibbs phenomenon, can also be observed in the numerical solutions of problems associated with flow through porous media with disparate properties. In layered porous domains, conventional continuous finite element methods are not capable of capturing abrupt changes in material properties and result in overshoots and undershoots in the profiles of numerical solutions along the interface of layers where there are jump discontinuities. In order to eliminate such erroneous oscillations, one possible approach is discontinuous Galerkin (DG) methods. DG methods have been successfully employed by

[Hughes et al., 2006] for the case of Darcy equations. An extension of the proposed framework using discontinuous Galerkin method for double porosity/permeability model can be obtained using a method similar to the one proposed by [Hughes et al., 2006]. However, obtaining such an extension and comparison between the performance of continuous and discontinuous formulations for capturing abrupt changes in material properties are beyond the scope of this paper and will be addressed in a subsequent one.A common assumption in models of flow in porous media is that of steady-state conditions. However, many flows occurring in porous media such as aquifers and oil-bearing strata are transient or unsteady in nature. In this paper, we extend the proposed stabilized mixed formulation for the double porosity/permeability mathematical model to the transient case, and we illustrate this extension can accurately capture the transient flow characteristics.

Recently, it has been shown that some stabilized methods (which are primarily designed to suppress numerical instabilities) when applied to solve problems with physical instabilities, suppress both types of instabilities [Shabouei and Nakshatrala, ; Shabouei, 2016]. Therefore, a good test of the proposed stabilized mixed formulation for a coupled flow and transport problem involves a problem that exhibits a physical instability similar to the classical Saffman-Taylor instability [Saffman and Taylor, 1958]. Using numerical simulations we show that the proposed formulation suppresses only the spurious numerical instabilities while capturing the underlying physical instability.

The rest of this paper is organized as follows. After an outline of the governing equations of the double porosity/permeability model in Section 2, the corresponding stabilized mixed formulation is presented in Section 3 with a derivation provided in Appendix A. The theoretical convergence analysis for the proposed stabilized mixed formulation is presented in Section 4, followed by the numerical convergence behavior of the elements presented in Section 5 where patch tests in one- and three-dimensional spaces are described. The representative numerical results are used to showcase the performance of the proposed mixed formulation in Section 6. Section 7 provides the mechanics-based assessment of the numerical accuracy. The transient analysis and the capability of the computational framework for modeling coupled problems and capturing well-known physical instabilities in fluid mechanics are discussed in Sections 8 and 9. Finally, conclusions are drawn in Section 10.

Throughout this paper, repeated indices do not imply
summation. The terms *classical mixed formulation*
and *Galerkin formulation* are used interchangeably.

## 2. Governing Equations for Double Porosity/permeability

For convenience to the reader and for future referencing, we document the equations that govern the double porosity/permeability mathematical model considered in [Nakshatrala et al., 2016]. Let be a bounded domain, where “” denotes the number of spatial dimensions. The boundary of the domain is assumed to be piecewise smooth. Mathematically, , where the superposed bar denotes the set closure [Evans, 1998]. A spatial point is denoted by . The gradient and divergence operators with respect to are denoted by and , respectively. The unit outward normal to the boundary is denoted by .

The porous domain is assumed to consist of two dominant pore-networks, which will be referred to as the macro-pore and micro-pore networks and are, respectively, denoted by subscripts and . These pore-networks are connected with the possibility of mass exchange between them. The pressure field and the discharge (or Darcy) velocity in the macro-pore network are, respectively, denoted by and , and the corresponding ones in the micro-pore network are denoted by and . The governing equations under the double porosity/permeability model take the following form:

(2.1a) | ||||

(2.1b) | ||||

(2.1c) | ||||

(2.1d) | ||||

(2.1e) | ||||

(2.1f) | ||||

(2.1g) | ||||

(2.1h) |

where is the specific body force. The true density and coefficient of viscosity of the fluid are, respectively, denoted by and .

denotes the permeability tensor for macro-pore (

) and micro-pore () networks. denotes the part of the boundary on which the normal component of the velocity is prescribed in the macro-pore () and micro-pore () networks. Similarly, is that part of the boundary on which the pressure is prescribed in the macro-pore () and micro-pore () networks. and denote the prescribed pressures on and , respectively. and denote the prescribed normal components of the velocities on and , respectively. is the rate of volume exchange of the fluid between the two pore-networks per unit volume of the porous medium, and we model it as follows [Barenblatt et al., 1960]:(2.2) |

where is a dimensionless characteristic of the porous medium. In the rest of the paper, as is commonly done in the literature, will be simply referred to as the mass transfer. For mathematical well-posedness, we assume that

(2.3) |

## 3. A Stabilized Mixed Weak Formulation

In this section, we present the proposed stabilized mixed formulation for the double porosity/permeability model. A derivation of the proposed formulation is provided in Appendix A. The proposed formulation is built upon the stabilization ideas put forth in a pioneering paper by [Masud and Hughes, 2002]. The proposed formulation for double porosity/permeability model can be obtained by adding a stabilization term, similar to the one proposed by [Masud and Hughes, 2002] for the case of single-pore network Darcy equations, to each pore-network. The stabilization terms are based on the residual of the balance of linear momentum in each pore-network. The stability can be achieved without adding residual-based stabilization terms related to the mass balance equations for any of the pore-networks. We also present an extension of the proposed formulation for enforcing the velocity boundary conditions weakly, which will be convenient for problems involving curved boundaries. This extension is achieved by employing a procedure similar to the one proposed by [Nitsche, 1971].

We define the relevant function spaces, which will be used in the rest of this paper. We denote the set of all square-integrable functions on by . For mathematical well-posedness, we assume that

(3.1) |

where is a non-integer Sobolev space [Adams and Fournier, 2003]. The function spaces for the velocity and pressures fields are defined as follows:

(3.2a) | ||||

(3.2b) | ||||

(3.2c) | ||||

(3.2d) | ||||

(3.2e) | ||||

(3.2f) |

where is a standard Sobolev space, and is the dual space corresponding to [Adams and Fournier, 2003]. The standard inner-product over a set is denoted as

(3.3) |

For convenience, the subscript will be dropped if

. Moreover, the action of a linear functional on a vector from its associated vector space is denoted by

.A few remarks are needed regarding the following condition on the pressures in the function spaces and :

This condition of vanishing mean pressure
in one of pore-networks is a mathematically
elegant way of fixing the datum for the
pressure. Without fixing the datum for
the pressure (which will be the case
when only the velocity boundary conditions
are prescribed on the entire boundary), one
can find the pressures only up to an arbitrary
constant, which will be the case even under Darcy
equations
[Nakshatrala et al., 2006].
Herein, we introduced the vanishing mean
pressure condition into the function spaces
to ensure uniqueness of the solutions, which
will be established later in this paper.
However, it should be emphasized that
vanishing mean pressure in one of the
pore-networks is *not* necessary
for all the problems under the double
porosity/permeability model.
One can fix the datum for the pressure
under the double porosity/permeability
model by prescribing the pressure in
at least one of the pore-networks on
a portion of the boundary, which is a
set of non-zero measure. To put it
differently, for problems with
pressure boundary conditions, the
datum for the pressure is automatically
fixed through the prescribed boundary
condition, and hence, for those problems,
one does not include the zero mean
pressure condition in the function
spaces and .
For example, see the problem in
subsection 5.1,
which deals with prescribed pressure boundary
conditions.

The *classical mixed formulation*, which is
based on the Galerkin formalism, reads as follows: Find
and such that we have

(3.4) |

where the bilinear form and the linear functional are, respectively, defined as follows:

(3.5) | ||||

(3.6) |

In a subsequent section, we will show that the equal-order interpolation for all the variables, which is computationally the most convenient, is not stable under the classical mixed formulation. Of course, one could use divergence-free elements (e.g., Raviart-Thomas spaces [Raviart and Thomas, 1977]) but they need special data structures and computer implementations. We, therefore, present a stabilized mixed formulation, which is stable under the equal-order interpolation for all the field variables.

[breakable]
The *proposed stabilized mixed formulation*
reads as follows: Find and such that we have

(3.7) |

where the bilinear form and the linear functional are, respectively, defined as follows:

(3.8) |

(3.9) |

In subsequent sections, we show that the proposed stabilized mixed formulation is consistent, stable and accurate.

### 3.1. Weak enforcement of velocity boundary conditions

In the previous derivations made earlier in this section, the pressure boundary conditions (i.e., equations (2.1g) and (2.1h)) are enforced weakly under the proposed stabilized mixed formulation and the classical mixed formulation. However, the velocity boundary conditions in which the normal components of the velocities are prescribed (i.e., equations (2.1e) and (2.1f)) are enforced strongly. For domains with curved boundaries, which are commonly encountered in subsurface modeling, it is desirable to even prescribe the velocity boundary conditions weakly. We, therefore, provide a possible extension of the proposed stabilized mixed formulation for weak enforcement of the velocity boundary conditions. To this end, we follow the approach proposed by [Nitsche, 1971]. The Nitsche’s method is a powerful tool for weakly enforcing Dirichlet boundary conditions without the use of Lagrange multipliers, and has been utilized by several works such as [Bazilevs and Hughes, 2007; Embar et al., 2010; Annavarapu et al., 2014; Schillinger et al., 2016]. The Nitsche’s method is sometimes referred to as a variationally consistent penalty method to enforce Dirichlet boundary conditions [Hansbo, 2005]. We extend the Nitsche’s method to the proposed four-field stabilized formulation to enforce the prescribed normal components of the velocities in the macro- and micro-pore networks.

The stabilized mixed formulation that enforces the velocity boundary conditions weakly can be obtained as follows: Find and such that we have

(3.10) |

where the bilinear form and the linear functional are, respectively, defined as follows:

(3.11a) | ||||

(3.11b) |

where is the mesh size and is the penalty parameter. In this paper, we have taken to be the maximum edge length in the mesh, and have taken the penalty parameter to be 10. In the above statement of the weak formulation, since the velocity boundary conditions are enforced weakly, the appropriate function space for the velocities and the associated weighting functions will be , which can be mathematically defined as follows:

(3.12) |

The function spaces for the pressures and their weighting functions, however, remain same as before (i.e., the space).

## 4. A Theoretical Analysis of the Proposed Mixed Formulation

In this section, we present a systematic mathematical analysis (i.e., existence, uniqueness and well-posedness) and error analysis (i.e., consistency, stability, order of convergence) of the proposed stabilized mixed formulation. For convenience, we define the following product spaces:

(4.1) |

We group the field variables as follows:

(4.2a) | ||||

(4.2b) |

Then, the proposed mixed formulation in equation (3.7) can be compactly written as: Find such that we have

(4.3) |

We shall establish the stability of the formulation under the following norm:

(4.4) |

where denotes the norm corresponding to the standard inner-product. We need to first show that is in fact a norm on and . To this end, the following lemma will be used.

###### Lemma 4.1.

(A property of semi-norms) If and are semi-norms, then is also a semi-norm.

###### Proof.

The homogeneity of directly stems from the homogeneity of the semi-norms and . To wit,

(4.5) |

The non-negativity of is straightforward; that is, . The triangle inequality for the semi-norms and implies that

(4.6) |

These inequalities imply that

(4.7) |

We have employed the AM-GM inequality in obtaining equation (4.7), which further implies that

(4.8) |

This establishes the triangle inequality for . The homogeneity, non-negativity and triangle inequality imply that is a semi-norm. ∎

###### Proposition 4.1.

(Stability norm) is a norm on and .

###### Proof.

We first note that and are symmetric and positive definite tensors. The square root of a symmetric and positive definite tensor exists, and is itself a symmetric and positive definite tensor [Gurtin, 1981]. This implies that the following individual terms form semi-norms on and :

(4.9) |

Then, Lemma 4.1 implies that is a semi-norm. It is easy to show that implies that

(4.10) |

where is a constant. Noting that and utilizing the following condition in the definition of :

(4.11) |

we conclude that . With this, we have established that = 0 implies that . Hence, is a norm. ∎

###### Theorem 4.1.

(Uniqueness of weak solutions) The weak solution under the proposed mixed formulation is unique.

###### Proof.

On the contrary, assume that and are both (weak) solutions of the weak formulation. This implies that

(4.12) |

By subtracting the above two equations and noting the linearity in the second slot, we obtain

(4.13) |

Since , we can choose . This particular choice implies that

(4.14) |

Using Proposition 4.1 (which establishes that is a norm on ) we conclude that . ∎

###### Theorem 4.2.

(Boundedness) The bilinear form is bounded. That is,

(4.15) |

where is a constant.

###### Proof.

A direct application of the triangle inequality of the absolute value on real numbers implies that

(4.16) |

Cauchy-Schwartz inequality on inner-product implies that

(4.17) |

By applying Cauchy-Schwartz inequality on -tuple real numbers (i.e., on Euclidean spaces) we obtained the following:

(4.18) |

That is, we have established that

(4.19) |

which completes the proof. ∎

###### Theorem 4.3.

(Coercivity) The bilinear form is coercive. That is, the bilinear form is bounded below.

###### Proof.

The coercivity of the bilinear form can be established from the definition of and Proposition 4.1 (i.e., is a norm on ) as

(4.20) |

∎

Given the coercivity and boundedness of the bilinear form and the continuity of the linear functional, one can conclude that the proposed mixed weak formulation is well-posed by invoking the Lax-Milgram theorem [Brenner and Scott, 1994].

### 4.1. Convergence and error analysis of the finite element formulation

We decompose the computational domain into “” subdomains (which will be the elements in the context of the finite element method) such that

(4.21) |

where a superposed bar indicates the set closure. We denote the finite element solution by . That is,

(4.22) |

Likewise,

(4.23) |

If we denote the set of all polynomials up to and including -th order over a set by , and the set of all continuous functions defined on (which is the set closure of ) by , then the following finite-dimensional spaces can be defined:

(4.24a) | ||||

(4.24b) | ||||

(4.24c) | ||||

(4.24d) | ||||

(4.24e) |

We define the corresponding product spaces as follows:

(4.25) |

It is important to note that and are closed linear subspaces of and , respectively. The finite element formulation corresponding to the proposed stabilized mixed formulation reads: Find such that we have

(4.26) |

In a given coordinate system, we denote . For a given multi-index (i.e., tuple) of non-negative integers, , with order , the corresponding partial derivative of a scalar field, , can be written as follows:

(4.27) |

Using the above notation, the Sobolev semi-norm, , for scalar and vector fields can be compactly written as follows:

(4.28) | |||

(4.29) |

where denotes the summation over all the possible tuples of non-negative integers with order , and denotes the characteristic length of the domain.

###### Remark 4.1.

Although the notation introduced in equation (4.27

) is common in the theory of partial differential equations (e.g.,

[Evans, 1998]), it may not be that common in the engineering literature. For the benefit of the reader, we provide the following few examples to make the notation more apparent: