{'fraction_non_alphanumeric': 0.03780237183958974, 'fraction_numerical': 0.012519089726427722, 'mean_word_length': 4.778407233903475, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On mechanics-based a posteriori criteria to assess accuracy of numerical solutions for Darcy and Darcy-Brinkman equations

In this paper, we present several mathematical properties that the solutions to Darcy and Darcy-Brinkman equations satisfy. These properties can serve as robust a posteriori error estimation techniques to verify numerical solutions for these equations. The mathematical properties include the total minimum mechanical power, the minimum dissipation theorem, a reciprocal relation, and a maximum principle for the vorticity. In particular, we show that, for a given set of boundary conditions, Darcy velocity has the minimum total mechanical power of all the kinematically admissible vector fields. We will also show that a similar result holds for a Darcy-Brinkman velocity. We also show that for a conservative body force, the Darcy velocity and the Darcy-Brinkman velocity have the minimum total dissipation among their respective kinematically admissible vector fields. Using numerical examples, we show that the minimum dissipation theorem can be utilized identify pollution errors in numerical solutions. We then show that the solutions to Darcy and Darcy-Brinkman equations satisfy a reciprocal relation, which has the potential to identify errors in the numerical implementation of boundary conditions. We also show that the vorticity under both steady and transient Darcy-Brinkman equations satisfy maximum principles if the body force is conservative body force and the permeability is homogeneous and isotropic. A discussion on the nature of vorticity under steady and transient Darcy equations is also presented. Using several numerical examples, we will demonstrate the predictive capabilities of the proposed a posteriori techniques in assessing the accuracy of numerical solutions for a general class of problems, which could involve complex domains and general computational grids.


INTRODUCTION
Porous media models are widely employed to make predictions in subsurface modeling, which typically involves complex geometries and heterogeneous material properties. It therefore becomes necessary in subsurface modeling to rely on numerical simulations to solve realistic problems so as to make critical engineering decisions. For a problem with no known analytical or reference solution (which is typically the case in subsurface modeling), it becomes a non-trivial task for a modeler to ensure that the employed numerical formulations and simulators are sufficiently accurate. Below are some specific challenges that a subsurface modeler may face in using numerical simulators: (i) How much mesh refinement is required to obtain solutions with desired degree of accuracy for a problem that does not have an analytical or reference solution?
(ii) Will the chosen mesh be able to resolve singularities in the solution and avoid pollution errors? To put it differently, can we identify whether a particular type of mesh suffers from pollution errors for a problem with singularities? (iii) Has the computer implementation been done properly? (iv) Is the chosen numerical formulation accurate/appropriate for the chosen problem?
In the literature, one finds the usual approach of employing tools from functional analysis to obtain a priori estimates and assess stability. For example, see [Brezzi and Fortin, 1991;Babuška and Strouboulis, 2001]. On the other hand, this paper aims to address the aforementioned questions for popular porous media models by providing various a posteriori techniques with firm mechanics underpinning. It needs to be emphasized that the techniques presented in this paper will not be able to completely address the aforementioned challenges. For example, there is currently no methodology that can judge with certainty the correctness of a computer implementation. The proposed techniques are no exception. However, the proposed a posteriori techniques can be judiciously utilized to increase the confidence or to detect shortcomings of a computer implementation, and assess the accuracy of numerical solutions.
1.1. Popular porous media models. The simplest and yet the most popular model that subsurface modelers use is the Darcy model [Darcy, 1856], which describes the flow of an incompressible fluid in a rigid porous media. The Darcy model has been employed in several branches of engineering and in several important technological applications. Some specific examples include ground-water hydrology [Raats, 1984;Munaf et al., 1993], enhanced oil recovery [Minkoff et al., 2003], and simulation of flow of resin in composite manufacturing Lee et al., 2009]. The Darcy model has also been a cornerstone for the theoretical development of mixture theories [Bowen, 1976;de Boer, 2000]. Darcy equations have posed challenges to the finite element community, and played a crucial role in the development of mixed and stabilized formulations [Brezzi and Fortin, 1991;Masud and Hughes, 2002;]. These equations have received considerable attention even from the mathematical community, especially in the area of mathematical homogenization [Hornung, 1996]. Brinkman [Brinkman, 1947a,b] proposed a popular extension to the Darcy model, which is commonly referred to as the Darcy-Brinkman model. In addition to drag between the fluid and the porous solid, the Darcy-Brinkman model accounts for the friction between fluid layers. The Darcy and Darcy-Brinkman models are adequate for most of the civil engineering and geotechnical applications in which the velocities and gradient of velocities are small. It is, however, not possible to obtain analytical solutions under these models, and one commonly seeks numerical solutions for realistic problems. Some notable numerical formulations for Darcy and Darcy-Brinkman equations are [Chavent et al., 1984;Ewing and Heinemann, 1984;Durlofsky, 1994;Arbogast et al., 1997;Masud and Hughes, 2002;Brezzi et al., 2005;Burman and Hansbo, 2007;Nakshatrala and Rajagopal, 2011;Srinivasan and Nakshatrala, 2012]. The question that still remains is how to verify the accuracy of numerical solutions, especially for realistic problems. This paper addresses this issue by providing several a priori techniques to verify the accuracy of numerical solutions, which falls under a current active area of research of validation and verification. (V&V). Errors can arise in both physical modeling and numerical simulation. The study of errors due to physical modeling is referred to as validation, and the study of error in a numerical simulation is referred to as verification. As Blottner [Blottner, 1990] nicely puts it, validation is to solve "right governing equations" and verification is to solve "governing equation right".

Validation and verification
Validation errors arise when a model is used out of its application range, and the boundary conditions are themselves an important source of validation errors. The validation is beyond the scope of this study, and the reader can consult [Roache, 1998] and the references therein. The errors in the verification can arise from three broad sources including numerical errors, round-off errors (due to the finite precision arithmetic), and programming mistakes [Oberkampf et al., 2004]. Basically, the verification is to ensure that the code produces a solution to the model with some degree of accuracy, and the numerical solution must be consistent. Verification itself can be divided into two categories: verification of code and verification of calculation [Roache, 1998].
1.2.1. Verification of code. Verification of code addresses the question of whether the numerical algorithms have been programmed and implemented correctly in the code. The two currently popular approaches to verify a code are the method of exact solutions (MES), and the method of manufactured solutions (MMS). In the method of exact solutions, one first seeks exact solutions often with simplifications to the equations and/or the boundary conditions. The numerical results with corresponding inputs are then compared to the exact solutions. In the method of manufactured solutions, an analytical solution is manufactured a priori, and a source term is then calculated that satisfies the governing equations. Steinberg and Roache [Steinberg and Roache, 1985] advocate the combined use of the method of manufactured solutions and grid convergence for the verification of code. If the code passes a desired number of test problems, which have been constructed by MES and MMS, then the code is certified to be verified adequately. This by no means implies that the code is completely free from bugs. A more thorough discussion on MES and MMS can be found in [Knupp and Salari, 2003;Roy et al., 2004;Roache, 1998].
1.2.2. Verification of calculation. Verification of calculation (which is also referred to as solution verification) is to estimate the overall magnitude (not just the order) of the numerical errors in a calculation, and the procedure invariably involves a posteriori error estimation [Salari and Knupp, 2000]. The numerical errors in the solution verification can arise from two different sources including discretization errors and solution errors. The discretization errors refer to all the errors caused by conversion of the governing equations (PDEs and boundary conditions) into discrete algebraic equations, whereas the solution errors refer to the errors in approximate solution of the discrete equations. The numerical errors include errors from insufficient mesh resolution, an improper selection of time-step, and an incomplete iterative convergence. For more details on verification of calculation, see [Roy et al., 2004;Salari and Knupp, 2000;Roache, 1997Roache, , 1998Babuška and Oden, 2004;Oberkampf and Blottner, 1998;Oberkampf et al., 2004].
1.3. A posteriori techniques. The aim of a posteriori error estimation is to assess the accuracy of the numerical approximation in the terms of known quantities, which include geometrical properties of computational grid, the input data, and the numerical solution. A posteriori error techniques monitor various forms of the error in the numerical solution such as velocity, stress, mean fluxes, and drag and lift coefficients [Becker and Rannacher, 2001]. Such error estimation differ from a priori error estimates in that the error controlling parameters depend on the unknown quantities. A priori error estimation investigates the stability and convergence of a solver and can give rough information on the asymptotic behavior of errors in calculations when grid parameters are changed appropriately [Ainsworth and Oden, 1997].
In the current study, we propose new mechanics-based criteria that can be used to verify accuracy and convergence of numerical solutions to Darcy and Darcy Brinkman equations. The solution verification requires confirmation of grid convergence, which is one of the most common and reliable error estimation methods [Roache, 1997]. Similar to other grid convergence studies, we address only solution verification. Since the aim of this paper is a posteriori error estimation, we assume that the code itself has already been verified for the Darcy and the Darcy-Brinkman class of problems so that the programming mistakes are not an issue. We are not also concerned with validation. It means that the Darcy and Darcy-Brinkman models are physically adequate to model the problems.
1.4. Organization of the paper. Section 2 presents the governing equations arising from the Darcy and Darcy-Brinkman models. In Section 3, we propose various mathematical properties that the solutions to these governing equations satisfy. We will also discuss how these properties can be utilized as robust a posteriori criteria to assess the accuracy of numerical solutions. Section 4 presents several steady-state numerical results to illustrate the predictive capabilities of the proposed a posteriori criteria with respect to singularities, pollution errors, and discretization errors in the implementation of (Neumann) boundary conditions. In Section 5, we shall use synthetic reservoir data to demonstrate the usefulness of the proposed a posteriori techniques, especially for problems involving spatially heterogeneous medium properties. Section 6 discusses a posteriori criteria for transient problems, and presents representative numerical results in support of the theoretical predictions. Finally, conclusions are drawn in Section 7.

DARCY AND DARCY-BRINKMAN MODELS
Let Ω ⊂ R nd be an open and bounded domain, where "nd" denotes the number of spatial dimensions. We shall denote the set closure of Ω by Ω. Let ∂Ω := Ω − Ω denote the boundary, which is assumed to be piecewise smooth. A spatial point in Ω is denoted by x. The spatial gradient and divergence operators are, respectively, denoted as grad [·] and div [·]. Let v : Ω → R nd denote the velocity field and p : Ω → R denote the pressure field. The symmetric part of the gradient of velocity is denoted by D(x). That is, The unit outward normal to the boundary is denoted as n(x). The boundary is divided into two parts: Γ v and Γ t . Γ v is the part of the boundary on which the velocity is prescribed, and Γ t is that part of the boundary on which the traction is prescribed. For mathematical well-posedness, The porous media models that will be considered in this paper are the Darcy and Darcy-Brinkman models. Both these models describe the flow of an incompressible fluid through rigid porous media. We completely neglect the motion of the porous solid. The Cauchy stress in the Darcy and Darcy-Brinkman models, respectively, take the following form: where I denotes the second-order identity tensor, and µ is the dynamic coefficient of viscosity. The steady-state governing equations based on the Darcy model can be written as follows: where α(x) is the drag coefficient, ρ is the density of the fluid, b(x) is the specific body force, v n (x) is the prescribed normal component of the velocity, and p 0 (x) is the prescribed pressure. The steady-state governing equations based on the Darcy-Brinkman model take the following form: where v p (x) is the prescribed velocity vector, and t p (x) is the prescribed traction. We shall call a vector-field to be Darcy velocity if satisfies equations (2.3a)-(2.3d). We shall call a vector field to be Darcy-Brinkman velocity if it satisfies equations (2.4a)-(2.4d). Equations (2.3a) and (2.4a) can be obtained from the balance of linear momentum under the mathematical framework offered by the theory of interacting continua [Nakshatrala and Rajagopal, 2011]. The drag term, α(x)v(x), models the frictional force between the fluid and the porous solid. The term div[2µD] in the Darcy-Brinkman model arises due to the internal friction between the layers of the fluid. The pressure, p(x), is an undetermined multiplier that arises due to the enforcement of the incompressibility constraint given by equations (2.3b) and (2.4b). The drag coefficient is related to the coefficient of viscosity of the fluid and the permeability, k(x), as follows: In general, it is not possible to obtain analytical solutions to the systems of equations given by either (2.3a)-(2.3d) or (2.4a)-(2.4d). Hence, one needs to resort to numerical solutions. This paper does not concern with developing new numerical formulations to solve the aforementioned mathematical models. The paper instead focuses on deriving mathematical properties with firm mechanics underpinning that the solutions to these mathematical models satisfy. We shall then illustrate how these mathematical properties can serve as robust "a posteriori" criteria to assess the accuracy of numerical solutions.

MATHEMATICAL PROPERTIES: STATEMENTS AND DERIVATIONS
In the remainder of this paper, we shall refer to a vector field v : Ω → R nd as kinematically admissible if it satisfies the following conditions: in Ω), and (ii) v(x) satisfies the boundary conditions.
It needs to emphasized that a kinematically admissible vector field need not satisfy the balance of linear momentum given by equation (2.3a) or (2.4a). Clearly, the Darcy velocity and the Darcy-Brinkman velocity are kinematically admissible vector fields. For some of the results presented in this paper, we will need the body force to be conservative, which is a common terminology in potential theory [Kellogg, 2010]. The body force ρb(x) is said to be conservative if there exists a scalar potential ψ(x) such that ρb(x) = −grad [ψ]. We shall define the dissipation functional as follows: Since α > 0 and µ > 0, it is straightforward to check that Φ[v] is a norm. In fact, it can be shown that Φ [v] under the Darcy model is equivalent to the natural norm in (L 2 (Ω)) nd , where (L 2 (Ω)) nd is a space of square integrable vector fields defined from Ω to R nd . Similarly, it can be shown that Φ [v] under the Darcy-Brinkman model is equivalent to the natural norm in (H 1 (Ω)) nd , which is a Sobolev space. For further details on function spaces and norms, refer to [Evans, 1998].
In this section, we shall present four important mathematical properties that the solutions to Darcy equations and Darcy-Brinkman equations satisfy. These properties will be referred to as (i) the minimum total mechanical power theorem, (ii) the minimum dissipation theorem, (iii) a reciprocal relation, and (iv) a maximum principle for vorticity. As a passing comment, we will employ the minimum dissipation theorem to show the uniqueness of solution for Darcy equations and Darcy-Brinkman equations. In the porous media literature, these results have been neither discussed nor utilized to solve problems. More importantly, these results have not been used to assess the accuracy and convergence of numerical solutions of porous media models. For example, it will be shown that the minimum total mechanical power theorem can be utilized to assess the accuracy of the implementation of both Dirichlet and Neumann boundary conditions. On the other hand, the minimum dissipation theorem can be used to identify pollution error in the numerical solutions. Herein, we give detailed mathematical proofs for Darcy-Brinkman equations. We however provide comments on the corresponding proofs for Darcy equations.
Theorem 1 (Minimum total mechanical power theorem). Let v(x) be the Darcy-Brinkman velocity vector field. Then, any kinematically admissible vector field v(x) satisfies the following inequality: That is, for given boundary conditions, body force and tractions; the Darcy-Brinkman velocity will have the minimum total mechanical power among all the possible kinematically admissible vector fields.
From the hypothesis of the theorem, δv(x) satisfies the following relations: δv(x) = 0 ∀x ∈ ∂Ω (3.5a) div[δv] = 0 ∀x ∈ Ω (3.5b) Let us start with the dissipation due to the vector field v(x): Using equations (2.4a) and (2.2b) the first integral in the above equation can be written as follows: The symmetry of D(x) allows the second integral to be written as follows: Noting that div[δv] = 0 in Ω and by employing Green's identity, we have the following: From equations (3.7) and (3.8), inequality (3.6) can be written as follows: This completes the proof.
Remark 1. For Darcy equations (2.3a)-(2.3d), one can state the minimum total mechanical power theorem as follows: For given boundary conditions, and body force; the Darcy velocity, v(x), has the minimum total mechanical power among all the kinematically admissible vector fields. That is, Theorem 2 (Minimum dissipation inequality). Let v(x) be the Darcy-Brinkman velocity vector field, and let Γ v = ∂Ω. Then, any kinematically admissible vector field v(x) satisfies the following inequality: (3.11) That is, for given velocity boundary conditions and conservative body force, the Darcy-Brinkman velocity has the minimum total dissipation due to drag and internal friction of all the possible kinematically admissible vector fields.
Proof. We shall employ the notation introduced in equation (3.4). Recall that the mechanical dissipation under the Darcy-Brinkman model is Let us start with the difference in total dissipation, and from inequality (3.6) we have: Using Green's identity, the first integral can be simplified as follows: Noting the symmetry of D(x) and using Green's identity, the total dissipation due to internal friction can be simplified as follows: From equations (3.7)-(3.8), the total dissipation due to drag and internal friction satisfies: This completes the proof.
The solutions to Darcy and Darcy-Brinkman equations posses reciprocal relations similar to the famous Betti's reciprocal relation in the theory of linear elasticity [Truesdell and Noll, 2004;Sadd, 2009] and to a classical reciprocal relation in the area of creeping flows [Guazzelli and Morris, 2012]. The Betti's reciprocal relation is often employed to solve a class of problems in linear elasticity, which otherwise may be difficult to solve. Mathematically, the Betti's reciprocal relation is equivalent to the existence and symmetry of Green's function. We now precisely state a reciprocal relation that the solutions of Darcy-Brinkman equations satisfy, and then provide a mathematical proof.
2 (x)}, respectively. Then, these fields satisfy the following relation: Proof. Let us start with the left side of equation (3.17). Noting that v p (x) = 0 on Γ v and Γ v ∪ Γ t = ∂Ω, one can proceed as follows: In the above step, we have used the fact that D 1 (x) is a symmetric second-order tensor. Since div [v 2 ] = 0 we have the following: Similarly, it can be shown that the right side of equation (3.17) is also equal to This completes the proof.
The following notation will be used later to verify the reciprocal relation: where the left and right integrals are, respectively, defined as follows: Remark 2. The corresponding reciprocal relation for Darcy equations can be written as follows: respectively. Then, these fields satisfy the following relation: Remark 3. It needs to be emphasized that the reciprocal relation will not directly be able to assess the accuracy of the pressure field in the computational domain. The reciprocal relation is ideal for assessing the accuracy of the velocity vector field in the domain and the accuracy of the implementation of prescribed traction boundary conditions. This reciprocal relation will not be able to provide information about the accuracy of the implementation of non-zero velocity boundary conditions.
Next, we discuss in the form of a theorem on the nature of the vorticity under the Darcy and Darcy-Brinkman models. To this end, In a Cartesian coordinate system, the components of vorticity take the following form: It should be emphasized that curl [·] operator is defined only in R 3 (i.e., the three-dimensional Euclidean space). However, for two-dimensional problems, one can consider the vorticity as follows: where z denotes the axis perpendicular to the two-dimensional plane in which the problem is defined, and e z is the unit vector along the z-direction.
Theorem 4 (On the nature of vorticity under Darcy and Darcy-Brinkman models). Assume that the medium is isotropic and homogeneous (i.e., α(x) is a constant scalar), the body force is a conservative vector field (i.e., ρb(x) = −grad[ψ(x)]), and the response is steady-state. Then the vorticity vanishes under Darcy equations. Under Darcy-Brinkman equations, the vorticity is an eigenvector of the Laplacian with 1/k as the eigenvalue. Moreover, the vorticity satisfies a maximum principle in which the non-negative maximum and the non-positive minimum occur on the boundary.
Proof. By taking the curl on both sides of the balance of linear momentum under Darcy equations (2.3a) we obtain the following: Since α is spatially homogeneous scalar, one can conclude that the vorticity vanishes under the Darcy model. The incompressibility constraint implies that the balance of linear momentum under the Darcy-Brinkman model can be written as follows: (3.25) By taking curl on both sides of the above equation, we get the following: where ∆ denotes the Laplacian operator and k is the permeability. The above equation is a vector eigenvalue problem in which the vorticity vector is the eigenvector, and 1/k is the corresponding eigenvalue. For two-dimensional problems, we have the following: where ∆ 2D denotes the two-dimensional Laplacian operator. The above equation is a scalar eigenvalue problem in which ω z (x, y) is the eigenvector and 1/k is the corresponding eigenvalue.
Since 1/k > 0, this equation is commonly referred to as diffusion with decay, which is a linear self-adjoint elliptic partial differential equation. It is well-known that such a partial differential equation satisfies a maximum principle [Gilbarg and Trudinger, 2001]. Mathematically, the maximum principle for the vorticity under the Darcy-Brinkman model takes the following form: If w z (x) ∈ C 2 (Ω) ∩ C 0 (Ω), then the non-negative maximum and the non-positive minimum occur on the boundary. That is, where C 2 (Ω) denotes the set of twice differentiable functions defined on Ω, and C 0 (Ω) is the set of functions that are continuous to the boundary.
The maximum principle can be utilized to verify the accuracy of numerical solutions by plotting vorticity, and checking whether the non-negative maximum and the non-positive minimum of the vorticity occur on the boundary. For the Darcy model with isotropic and homogeneous medium properties, and conservative body force, it can be shown that the vorticity vanishes (i.e., ω(x) = 0). It should be, however, noted that heterogeneity, pressure-dependent viscosity, or nonconservative body force can introduce vorticity under the Darcy model. All the above results can serve as invaluable tools to assess the performance of a numerical formulation to verify a computer implementation, and serve as metrics for numerical convergence.

STEADY-STATE NUMERICAL RESULTS
We shall first non-dimensionalize the governing equations by choosing primary variables that seem appropriate for problems arising in modeling of flows through porous media. This nondimensional procedure is different from the standard non-dimensionalization procedure for incompressible Navier-Stokes in the choice of primary variables. In the standard non-dimensionalization of Navier-Stokes equations, one employs characteristic velocity v, characteristic length L and density of the fluid ρ as primary variables. We shall choose L (reference length in the problem), g (acceleration due to gravity) and p atm (atmospheric pressure) as the reference quantities. Using these reference quantities, we define the following non-dimensional quantities: where the non-dimensional quantities are denoted by superposed bars. The gradient and divergence operators with respect to x are denoted as grad[·] and div[·], respectively. The scaled domain Ω scaled is defined as follows: A point in space with position vector x ∈ Ω scaled corresponds to the same point with position vector given by x = xL ∈ Ω. Similarly, one can define the scaled boundaries: ∂Ω scaled , Γ v scaled , and Γ p scaled . Using the above non-dimensionalization procedure, Darcy-Brinkman Table 1. Body force problem: Non-dimensional parameters used in the problem.
Parameter Value α 1 µ 1 and 0.001 equations can be written as follows: Similarly, one could write the corresponding non-dimensional form of Darcy equations. For simplicity, the "over-lines" will be dropped in the remainder of the paper. We shall use several flow through porous media problems with diverse boundary conditions to illustrate that the proposed dissipation, total mechanical power, and reciprocal relation can be used as a good measure for the accuracy and convergence of numerical results. We also show that the vorticity can be used to verify the accuracy of solution for the Darcy-Brinkman model.
We have employed P2P1 (which is based on six-node triangle element T6) and Q2P1 (which is based on nine-node quadrilateral element Q9) interpolations available in COMSOL [COM, 2013]. Typical structured finite elements utilized in this paper are shown in Figure 1. Unless mentioned otherwise, all the elements in a quadrilateral mesh are squares and all the elements in a triangular mesh are right-angled isosceles triangles. In this numerical solution study, we have taken h (the maximum element size) to be equal to the length of the side for square elements, and to the length of the base (or height) for right-angled isosceles triangles.
4.1. Body force problem. The test problem is pictorially described in Figure 2. The nondimensional parameters used in the numerical simulation are provided in Table 1. The conservative body force is taken as ρb(x) = 10[sin(πx), cos(πy)] (i.e., a = 10). Velocity boundary condition is prescribed on the entire boundary (i.e., Γ v = ∂Ω). For the Darcy-Brinkman model, we assume v p (x) = 0; and for the Darcy model, we assume v n (x) = 0. Figure 3 shows the minimum total mechanical power and the minimum dissipation with mesh refinement. The numerical results for the reciprocal relation with mesh refinement are shown in Figure 4. All the numerical results are in accordance with the theoretical predictions for this test problem.
4.2. Lid-driven cavity problem. The two-dimensional lid-driven cavity problem is a benchmark study widely used to investigate the accuracy of numerical formulations for various fluid models [Burggraf, 1966;Erturk, 2009]. Figure 5 provides a pictorial description of the problem. The domain of the problem is a bi-unit square. Velocity boundary condition are prescribed on the entire boundary (i.e., Γ v = ∂Ω), which implies that the minimum dissipation theorem is applicable to this problem. It should be noted that lid-driven cavity problem is not compatible with Darcy equations, which need only the normal component of the velocity to be prescribed on the boundary. Table 2. Lid-driven cavity problem: Non-dimensional parameters used in the problem.
Parameter Value However, the lid-driven cavity problem demands the prescription of the entire velocity vector field on the boundary. We shall therefore employ Darcy-Brinkman model in our numerical simulations. It is crucial to note that the solution to the lid-driven cavity problem has singularities at the top corners, which arise due to velocity discontinuity on the boundary [Botella and Peyret, 1998;Batchelor, 2000]. Discontinuities and singularities commonly occur in fluid dynamics, solid mechanics, and structural dynamics; and it is a subject of intense research and grave concern in error estimation [Oberkampf et al., 2004]. Some representative works on pollution errors due to singularities and discontinuities are [Babuška and Oh, 1987;Babuška et al., 1995Babuška et al., , 1997Oden et al., 1998;Roache, 1998;Botella and Peyret, 2001]. But none of these studies utilized the minimum dissipation theorem to detect pollution errors in porous media models.
Herein, we use an adaptive mesh to resolve the singularities in the solution. The non-dimensional parameters used in the lid-driven cavity problem are presented in Table 2. Figures 6(a) and 6(b), respectively, show the uniform structured mesh and the adaptive mesh. For the adaptive mesh, we generate fine grid in the region close to top lid (see Figure 6(b)). Figures 6(c) and 6(d) show variation of the minimum dissipation with mesh refinement for these grids. Under the adaptive mesh, the dissipation decreased uniformly and converged to a constant value with mesh refinement. On the other hand, the dissipation increased monotonically with mesh refinement under the structured mesh. More importantly, due to the presence of singularities and pollution error, the dissipation did not hit a plateau even for very fine meshes (i.e., 1/h ≥ 220). The dissipation reached a plateau relatively quickly under the adaptive mesh (say 1/h = 10). This problem clearly illustrates that the minimum dissipation theorem can be used to identify pollution errors in numerical solutions, which is one of the main findings of this paper.
4.3. Pipe bend problem. As another application of the proposed techniques, we consider an engineering problem commonly found in the fluid mechanics literature which is called pipe bend problem (for example see References [Borrvall and Petersson, 2003;Hansen et al., 2005;Aage et al., 2008;Pingen et al., 2009;Challis and Guest, 2009;Hassine, 2012]). In the pipe bend problem the computational domain Ω is a square with side length L = 1. The classical problem has only the velocity boundary condition (i.e., Γ v = ∂Ω). However, to test our proposed solution verification methods, we employed our techniques not only to the velocity boundary condition (classical) problem (Figure 7), but also to velocity-pressure boundary condition problems ( Figures  9 and 11). 4.3.1. Velocity boundary condition. For the velocity boundary condition problem (Figure 7) we shall assume that Γ v = ∂Ω. An inflow parabolic velocity is enforced on a part of the left boundary and an outflow parabolic velocity on a segment of the bottom. Each the parabolic velocity profile has the unit maximum value (i.e., v xmax = 1 or v ymax = 1). Elsewhere, the homogeneous velocity is prescribed (i.e., v p (x) = 0 for the Darcy-Brinkman model and v n (x) = 0 for the Darcy model). The Table 3. Pipe bend problem: Non-dimensional parameters of the parabolic velocity boundary condition problem.

Parameter
Value α 1 and 10 µ 1 and 0.001 velocity boundary condition makes the problem compatible with the total mechanical power and the dissipation theorems. The non-dimensional parameters used in the problem are presented in Table 3. Figure 8 shows the variation of minimum total mechanical power and minimum dissipation with mesh refinement under the quadrilateral and triangular elements. The result of the numerical solutions verification are presented in Figure 8(a) and 8(b) for the minimum dissipation and the total mechanical power, respectively. The numerical error decreases and converges uniformly. 4.3.2. Parabolic velocity-pressure boundary condition. We shall assume a separate velocity and pressure boundary condition for this problem (Figure 9). The traction is prescribed on a part of the bottom boundary (i.e., t p (x) = −p atm n(x) on Γ t ). The velocity has parabolic profile with unit maximum value (v xmax = 1) prescribed on a segment of the left boundary. Elsewhere, the homogeneous velocity is enforced. On account of the traction and parabolic velocity boundary conditions, current problem is not compatible with the minimum dissipation and reciprocal theorems. It is only compatible with the total mechanical power theorem (see Theorems 2 and 3 assumptions). It should be noted that the solution to the problem has the singularity at near corners of the outlet (i.e., Γ t ). Hence, we again use an adaptive mesh to resolve the pollution in the solution. The non-dimensional parameters using in the problem are provided in Table 3. Figure 10 depicts the variation of the minimum total mechanical power with mesh refinement for the quadrilateral and triangular elements. The top figures show the uniform structured and adaptive meshes. Under the structured mesh, the total mechanical power increased uniformly with mesh refinement due to the polluted area in the computational domain. The error decreased uniformly and converged to a constant value with mesh refinement under the adaptive mesh. So, in addition to error estimation capability, the total mechanical power theorem can be used to identify the pollution errors in the numerical solutions. 4.3.3. Velocity-pressure boundary condition. The pictorial description of the problem is shown in Figure 11. The traction is prescribed on a segment of the left boundary (i.e., t p (x) = −p inj n(x) on Γ t 1 ) and bottom boundary (i.e., t p (x) = −p atm n(x) on Γ t 2 ). Elsewhere, the homogeneous velocity is enforced. The non-dimensional parameters used in the problem are presented in Table  4. Current problem is compatible for the reciprocal relation, due to the prescription of the traction and homogeneous velocity on the boundary. The variation of the error in the reciprocal relation with mesh refinement under the triangular and quadrilateral elements is shown in Figure 12. Due to the singularity near the corners of outlet (i.e., Γ t 2 ) the results (i.e., ε reciprocal ) under the regular structured mesh is uniform but increasing. However, the error decreases uniformly under the adaptive mesh.  Figure 13 provides a pictorial description of the problem. The non-dimensional parameters used in the numerical simulation are provided in Table 5. The domain is a W × L rectangle. The homogeneous velocity boundary condition is enforced on the top and bottom sides of the boundary. The traction is prescribed on the left side (i.e., t p (x) = −p inj n(x) on Γ t 1 ) and on the right side (i.e., t p (x) = −p atm n(x) on Γ t 2 ). We shall use this test problem to assess the accuracy of numerical solutions with respect to the reciprocal relation. Figure 14 depicts the variation of the error in the reciprocal relation with mesh refinement for the triangular and quadrilateral grids. The error in the reciprocal relations under Darcy equations is very close to zero for all the meshes. The error under Darcy-Brinkman equations decrease uniformly with mesh refinement. All the numerical results are in accordance with the theoretical predictions.

Pressure driven problem.
The domain is a square with side length L = 1. A pictorial description of the problem is given in Figure 15. The traction is prescribed on the left boundary (i.e., t p (x) = −p inj n(x) on Γ t 1 ) and on the middle of the right boundary (i.e., t p (x) = −p atm n(x) on Γ t 2 ). Elsewhere, homogeneous velocity boundary condition is enforced. Due to the prescription of the traction on a part of the boundary, this problem is not compatible with the minimum dissipation theorem, but is compatible with the reciprocal relation and the total mechanical power theorem. Herein, we present the results for the reciprocal relation. The non-dimensional parameters used in the numerical simulation are provided in Table 6. Figure 16 show the variation of the error in the reciprocal relation with mesh refinement for triangular and quadrilateral meshes. Since the problem has a singularity near the corners of Γ t 2 , the convergence under uniform structured meshes is slow which improved by employing an adaptive mesh. Table 6. Pressure driven problem: Non-dimensional parameters used in the problem.

Parameter
Value α 1 µ 1 and 0.001 ρ 1 p inj 5 and 7.5 p atm 1 L 1 4.6. Vorticity results. The maximum principle given by Theorem 4 can be used to assess the accuracy of numerical solutions to the Darcy-Brinkman equations by checking whether the non-negative maximum and the non-positive minimum of the vorticity occur on the boundary. Figure 17 shows that the maximum principle is satisfied for various problems under the steadystate Darcy-Brinkman equations.

SYNTHETIC RESERVOIR DATA: MARMOUSI DATASET
We will now solve an idealized reservoir problem using a popular synthetic dataset -the socalled (smooth) Marmousi dataset [Versteeg and Grau, 1990;Versteeg and Lailly, 1991;Versteeg, 1993;Klimes, 2014;Benamou, 2014]. The dataset provides spatially varying speed of sound on a 384 × 122 grid. We have assumed that the permeability scales linearly with the values provided by the dataset. This is just an arbitrary choice to generate a heterogeneous dataset for permeability. However, it should be noted that the conclusions that will be drawn here will be valid even if one uses another dataset for the permeability. Figure 18 shows the contours of Marmousi dataset.
The boundary value problem of the reservoir is pictorially described in Figure 19. This computational domain has been employed in some recent works (e.g., [Nakshatrala and Rajagopal, 2011]). All these works have assumed homogeneous medium properties, and did not use a reservoir data like the Marmousi dataset. Moreover, these studies did not address the use of a posteriori techniques to access numerical accuracy, which is the main focus of the current paper. The parameters used in this problem are provided in Table 7. Figure 20 shows that the errors in the reciprocal relation are larger under uniform structured meshes than under adaptive meshes. This is due to the fact that uniform structured meshes suffer from pollution errors due to the singularity near the production well. The variation of the error in the reciprocal relation with mesh refinement provides guidelines on how much refinement is required to obtain solution of some desired accuracy especially in those cases where there are no analytical solutions. Figure 21 shows the magnitude of the vorticity, and one can see from the figure that the maximum principle for the vorticity is satisfied.
in Ω (6.1e) Of course, the convective term grad[v]v is neglected. We assume that the coefficient of viscosity of the fluid and the permeability of the porous solid to be constants, and hence the drag coefficient is constant. We further assume that the density is homogeneous, and the body force is assumed to be conservative. The unsteady Darcy-Brinkman equations can be written as follows: where v p (x, t) is the prescribed velocity vector, and t p (x, t) is the prescribed traction.

Mathematical properties.
Under unsteady Darcy equations, the vorticity satisfies the following equation: which is an ordinary differential equation at each spatial point. The solution takes the following form: ω(x, t) = ω 0 (x) exp − αt ρ (6.4) where ∆ is the Laplacian operator. The above equation is a homogeneous linear parabolic partial differential equation, which is also known to satisfy a maximum principle. This implies that both the maximum and the minimum will occur either in the initial condition or on the boundary. One can check whether numerical solutions satisfy the aforementioned maximum principle. Also, whether refining the grid spacing and time-steps affect the performance of numerical solutions with respect to this metric. One can devise any test problem as long as the following assumptions are met: (i) µ is constant, (ii) permeability is homogeneous, and (iii) the body force is conservative.  Figure 13 provided a pictorial description of the pressure slab problem.) The non-dimensional parameters used in the numerical simulation are presented in Table 5, and the initial condition is provided in Table  8. Figure 22 shows that log( ω ω 0 ) is a straight line with slope equal to α ρ under the mesh and time refinements for various spatial points in the solution of the transient Darcy equations (for the current problem α ρ = −1). Figure 23 verifies the maximum principle for vorticity for various test problems under transient Darcy-Brinkman equations. In all the cases, the non-negative maximum and nonpositive minimum of the vorticity occur on the boundary, which agree with the mathematical theory presented above.

CONCLUDING REMARKS
We presented various test problems, which can serve as benchmark problems for verifying numerical implementation of solvers for the Darcy and Darcy-Brinkman models. We also presented four important properties (minimum dissipation theorem, minimum total mechanical power, reciprocal relation, and maximum principle for vorticity) that the solutions under the Darcy and the Darcy-Brinkman models satisfy. These properties can be effectively used to assess the accuracy of numerical solutions. An attractive feature is that these properties can be verified for any given problem (i.e., needed not be one of the benchmark problems), and for any computational domain.
For example, if the problem involves prescribing velocity boundary condition on the entire boundary, then one plots the dissipation with respect to mesh refinement. Some of the main conclusions are: (a) If the numerical formulation is converging, the dissipation, total mechanical power, or reciprocal relation should decrease with mesh refinement. If this does not occur, one needs to suspect that there are singularities in the solutions, or the numerical formulation does not perform well with respect to local mass balance property. (b) It has been shown that the minimum dissipation theorem can be utilized to identify pollution errors in numerical solutions. The theorem can also be used to assess whether a given type of mesh will be able to resolve singularities in the solution. This can be assessed by creating a series of hierarchical meshes, and plot the dissipation with respect to characteristic mesh size. A given type of mesh will resolve singularities in the solution and will not be affected by pollution errors if the total dissipation decreases uniformly and reaches a plateau with a hierarchical mesh refinement. (c) The non-negative maximum vorticity and the non-positive minimum vorticity under Darcy-Brinkman equations with homogeneous isotropic medium properties should occur on the boundary.
The proposed a posteriori techniques can be invaluable addition to the usual repertoire of methods for verification -the method of exact solutions and the method of manufactured solutions.

APPENDIX: UNIQUENESS OF SOLUTION
For completeness and pedagogical value, we now show that the uniqueness of solution under Darcy-Brinkman equations is a direct consequence of the minimum dissipation inequality. One can construct a proof for the uniqueness of solution under Darcy equations on similar lines.
Theorem 5 (Uniqueness theorem). The solution to Darcy-Brinkman equations (2.4a)-(2.4d) is unique up to an arbitrary constant for the pressure given b(x), v p (x) and t p (x).
Proof. On the contrary assume that {v 1 (x), p 1 (x)} and {v 2 (x), p 2 (x)} are two solutions to Darcy-Brinkman equations for the prescribed data. Let us consider the following quantity: Noting that div[v 1 ] = 0 and div[v 2 ] = 0, the second integral can be simplified as follows: In obtaining the above equation, we have used the fact that T 1 (x) and T 2 (x) are symmetric. Using Green's identity, the above equation can be written as follows: Using boundary conditions and the balance of linear momentum, we get the following: This implies that I = 0. Since α(x) > 0 ∀x ∈ Ω and µ > 0 in Ω, one can conclude that v 1 (x) = v 2 (x) ∀x ∈ Ω (7.1a) That is, the velocity and symmetric part of velocity vector field are unique. Using the equation for the balance of linear momentum (2.4a) and the fact that the velocity vector field is unique, one can obtain the following equation: This further implies that where p 0 is an arbitrary constant. This completes the proof.
Location of pressure unknown (c) Q2P1 interpolation for nine-node quadrilateral (Q9) element • Location of velocity unknown (d) P2P1 interpolation for six-node triangular (T6) element Figure 1. This figure shows the typical structured finite element meshes employed in the current study. We use Q2P1 and P2P1 mixed interpolations to approximate the unknowns (i.e., second-order interpolation for the velocity field, and first-order interpolation for the pressure field).   Table 1. The figure show that the total mechanical power and the dissipation decrease uniformly with mesh refinement, which are in accordance with Theorems 1 and 2.  x Figure 5. Lid-driven cavity problem: The computational domain is a unit square. Velocity boundary condition is prescribed on the entire boundary (i.e., Γ v = ∂Ω). The prescribed velocity on the top side is v x = 1 and v y = 0, and the prescribed velocity on the remaining sides of the boundary is zero (i.e., v p (x) = 0). Note that this test problem is not compatible with Darcy equations.
Results for Lid-driven cavity problem (a) Uniform structured mesh with 1/h = 12.   Table 2. One should note that the solution exhibits singularities at the top corners. It is evident that the convergence is uniform under both structured and adaptive meshes. However, the dissipation under the structured mesh increased with mesh refinements, whereas the trend is reversed under the adaptive mesh. This is because the adaptive mesh resolves the singularities, and removes the associated pollution error. This figure corroborates one of the main conclusions of the paper that the minimum dissipation theorem can be utilized to identify pollution errors in numerical solutions.    Table 3. Note that the total mechanical power is defined in equation (3.3). Elsewhere, the velocity is assumed to be zero (i.e., v p (x) = 0 for Darcy-Brinkman equations and v n (x) = 0 for Darcy equations). Note that the corners at the outlet are re-entrant corners.

Results for Pipe bend problem
(a) Uniform structured mesh with 1/h = 15.
(b) Adaptive mesh near the outlet with 1/h = 15 elsewhere.   Table 3. The total mechanical power increased uniformly with mesh refinement under the structured mesh, while the total mechanical power decreased uniformly and reached a plateau with mesh refinement under the adaptive mesh. The numerical results clearly demonstrate that the total mechanical theorem can be utilized to identify pollution errors due to singularities, and to assess whether a particular type of computational mesh is suitable for a given problem. Figure 11. Pressure pipe bend problem (pressure-pressure boundary condition): A pictorial description of the problem. The computational domain Ω is a unit square. The traction is prescribed on a part on the left side (i.e., t p (x) = −p inj n(x) on Γ t 1 ) and on a part on the bottom side (i.e., t p (x) = −p atm n(x) on Γ t 2 ). Homogeneous velocity is enforced on the remaining parts of the boundary.   Table 4. The top figures are based on the uniform structured mesh, and the bottom figures are based on the adaptive mesh. The adaptive mesh used in this problem is similar to the one provided in Figure 10(b). Under adaptive meshes, the convergence is uniform and the error decreases with mesh refinement. However, the convergence is slow, which is because of the discretization error associated with Neumann boundary conditions. Note the traction boundary condition (i.e., the Neumann boundary condition) is enforced weakly. Figure 13. Pressure slab problem: A pictorial description of the problem. The computational domain is a W × L rectangle. The traction is prescribed on the left side of the boundary (i.e., t p (x) = −p inj n(x) on Γ t 1 ) and on the right side (i.e., t p (x) = −p atm n(x) on Γ t 2 ). Elsewhere, homogeneous velocity is enforced.  Table 5. One can see the error in the numerical results with respect to the reciprocal relation for Darcy equations is very close to zero for all meshes. The corresponding error under Darcy-Brinkman equations decreases uniformly with mesh refinement.
Results for Pressure slab problem Figure 15. Pressure-driven problem: A pictorial description of the problem. The computational domain Ω is a square with length of each side L = 1. The traction is prescribed on the left side of the boundary (i.e., t p (x) = −p inj n(x) on Γ t 1 ) and on a part of the right side (i.e., t p (x) = −p atm n(x) on Γ t 2 ). Elsewhere, homogeneous velocity is enforced.
Results for Pressure-driven problem   Table 6. Clearly, the adaptive mesh performs better than the uniform structured mesh, and the minimum dissipation theorem can utilized to obtain information on the numerical performance of computational meshes even for a problem that does not possess an analytical solution.  [Benamou, 2014]. The dataset provides values on a 384 × 122 grid, which we have scaled for our rectangular computational domain of L x = 2 and L y = 1. This spatially heterogeneous dataset is widely used as a benchmark dataset in reservoir modeling.

Production well Injection well
Injection well Figure 19. Synthetic reservoir problem: This figure provides a pictorial description of the test problem. The domain of the problem is a rectangle of size H × L.
The injection pressure is prescribed on the left and right boundaries (i.e., t p (x) = −p inj n(x) on Γ t 2 ), and the atmosphere pressure is prescribed on the middle of the top side (i.e., t p (x) = −p atm n(x) on Γ t 1 ). The prescribed velocity on the remaining parts of the boundary is zero (i.e., v p (x) = 0 for the Darcy-Brinkman model and v n (x) = 0 for the Darcy model).   Table 7. The figures clearly show that the errors in the reciprocal relation are larger under the uniform structured meshes. This can be attributed to pollution errors due to singularities at the corners of the production wells.  Table 7. This figure verifies the maximum principle for the vorticity for the synthetic reservoir problem using Marmousi dataset.  The results show that the slope of log( ω ω 0 ) for various spatial points in the computational domain are close to − α ρ = −1. The figure also shows the slope is follows the theoretical prediction for various mesh refinements h and time-steps dt. The non-dimensional parameters used in this numerical experiment are provided in Tables 5 and 8.