Geostatistical modelling of soil contamination with arsenic, cadmium, lead, and nickel: the Silesian voivodeship, Poland case study

1.   Introduction

For simplicity, we consider Poisson equation with a Dirichlet boundary condition as our model problem.

where $\Omega$ is a bounded polygonal domain in $\mathbb{R}^2$.

Using integration by parts, we can get the variational form: find $u \in H^1(\Omega)$ satisfying $u = g\; \text{on}\; \partial\Omega$ and

Various finite element methods have been introduced to solve the Poisson equations (1)-(2), such as the Galerkin finite element methods (FEMs)[2, 3], the mixed FEMs [15] and the finite volume methods (FVMs) [6], etc. The FVMs emphasis on the local conservation property and discretize equations by asking the solution satisfying the flux conservation on a dual mesh consisting of control volumes. The mixed FEMs is another category method that based on the variable $u$ and a flux variable usually written as $p$.

The classical conforming finite element method obtains numerical approximate results by constructing a finite-dimensional subspace of $H_0^1(\Omega)$. The finite element scheme has the same form with the variational form (3): find $u_h \in V_h\subset H^1(\Omega)$ satisfying $u_h = I_hg\; \text{on}\; \partial\Omega$ and

where $V_h^0$ is a subspace of $V_h$ that satisfying $v_h = 0$ on $\partial\Omega$ and $I_h$ is the kth order Lagrange interpolation operator. The FE method is a popular and easy-to-implement numerical scheme, however, it is less flexible in constructing elements and generating meshes. These limitations are mainly due to the strong continuity requirements of functions in $V_h$. One solution to circumvent these limitations is using discontinuous approximations. Since the 1970th, many new finite element methods with discontinuous approximations have been developed, including the early proposed DG methods [1], local discontinuous Galerkin (LDG) methods [8], interior penalty discontinuous Galerkin (IPDG) methods [9], and the recently developed hybridizable discontinuous Galerkin (HDG) methods [7], mimetic finite differences method [10], virtual element (VE) method [4], weak Galerkin (WG) method [19, 20] and references therein.

One obvious disadvantage of discontinuous finite element methods is their rather complex formulations which are often necessary to ensure connections of discontinuous solutions across element boundaries. For example, the IPDG methods add parameter depending interior penalty terms. Besides additional programming complexity, one often has difficulties in finding optimal values for the penalty parameters and corresponding efficient solvers. Most recently, Zhang and Ye [21] developed a discontinuous finite element method that has an ultra simple weak formulation on triangular/tetrahedal meshes. The corresponding numerical scheme can be written as: find $u_h \in \tilde{V}_h$ satisfying $u_h = I_hg\; \text{on}\; \partial\Omega$ and

where $\tilde{V}_h$ is the DG finite element space and $\nabla_w$ is the weak gradient operator. The notion of weak gradient was first introduced by Wang and Ye in the weak Galerkin (WG) methods [19, 20]. The WG methods allow the use of totally discontinuous functions and provides stable numerical schemes that are parameter-independent and free of locking [17] in some applications. Another key feature in the WG methods is it can be used for arbitrary polygonal meshes. The WG finite element method has been rapidly developed and applied to other problems, including the Stokes and Navier-Stokes equations [11, 18], the biharmonic [14, 13] and elasticity equations [12, 17], div-curl systems and the Maxwell's equations and parabolic problem [23], etc. The introduction of the weak gradient operator in the conforming DG methods makes the scheme (5) maintain the simple formulation of conforming finite element method while have the flexibility of using discontinuous approximations. Hence, the programming complexity of this conforming DG scheme is significantly reduced. Furthermore, the scheme results in a simple symmetric and positive definite system.

Following the work in [21, 22], we propose a new conforming DG finite element method on rectangular partitions in this work. It can be obtained from the conforming formulation simply by replacing $\nabla$ by $\nabla_w$ and enforcing the boundary condition strongly. The simplicity of the conforming DG formulation will ease the complexity for implementation of DG methods. We note that the conforming DG method in [21] is based on triangular/tetrahedal meshes. Then in [22], the method is extended to work on general polytopal meshes by raising the degree of polynomials used to compute weak gradient.

In this paper, we keep the same finite element space as DG method, replace the boundary function with the average of the inner function, and use the weak gradient arising from local Raviart-Thomas (RT) elements [5] to approximate the classic gradient. Moreover, the derivation process in this paper is based on rectangular RT elements [16]. Error estimates of optimal order are established for the corresponding conforming DG approximation in both a discrete $H^1$ norm and the $L^2$ norm. Numerical verifications have been performed on different kinds of quadrangle finite element space. In particular, super-convergence phenomenon have been observed for $Q_0$ elements.

The rest of this paper is organized as follows: In Section 2, we shall present the conforming DG finite element scheme for the Poisson equation on rectangular partitions. Section 3 is devoted to a discussion of the stability and solvability of the new method. In Section 4, we shall prepare ourselves for error estimates by deriving some identities. Error estimates of optimal order in $H^1$ and $L^2$ norm are established in Section 5. In Section 6, we present some numerical results to illustrate the theory derived in earlier sections. Finally in section 7, we conclude our major contributions in this article.

Throughout this paper, we adopt the standard definition of Sobolev space $H^s(\Omega)$. For any given open bounded domain $K \subseteq \Omega$, $(\cdot, \cdot)_{s,K}, \|\cdot\|_{s,K}$, and $|\cdot|_{s,K}$ are used to denote the inner product, norm and semi-norm, respectively. The space $H^0(K)$ coincides with $L^2(K)$, and the subscripts $K$ in the inner product, norm, and semi-norm can be dropped in the case of $K = \Omega$. In particular, the function space $H^1_0(\Omega)$ is defined as

and the space $H(div, \Omega)$ is defined as the set of vector-valued functions $\boldsymbol{q}$, which together with their divergence are square integrable, i.e.

2.   Conforming DG method

Assume that the domain $\Omega$ is of polygonal type and is partitioned into non-overlapping rectangles $\mathcal{T}_h = \{T\}$. For each $T\in \mathcal{T}_h$, denote by $T^0$ its interior and $\partial T$ its boundary. Denote by $\mathcal{E}_h = \{e\}$ the set of all edges in $\mathcal{T}_h$, and $\mathcal{E}_h^0 = \mathcal{E}_h\setminus \partial\Omega$ the set of all interior edges in $\mathcal{T}_h$. For each $T\in\mathcal{T}_h$ and $e\in\mathcal{E}_h$, denote by $h_T$ and $h_e$ the diameter of $T$ and $e$, respectively. $h = \max_{T\in \mathcal{T}_h}h_T$ is the meshsize of $\mathcal{T}_h$.

For any interior edge $e\in\mathcal{E}^0_h$, let $T_1$ and $T_2$ be two rectangles sharing $e$, we define the average $\{\cdot\}$ and the jump $[\![ \cdot]\! ]$ on $e$ for a scalar-valued function $v$ by

where $v|_{\partial T_{i}},i = 1,2$ is the trace of $v$ on $\partial T_i$, $\boldsymbol{n</italic>}_1$ and $\boldsymbol{n</italic>}_2$ are the two unit outward normal vectors on $e$, associated with $T_1$ and $T_2$, respectively. If $e$ is a boundary edge, we define

We define a discontinuous finite element space

and its subspace

where $Q_k(T), \;k \geq 1$ denotes the set of polynomials with regard to quadrilateral elements. The weak gradient for a scalar-valued function $v\in V_h$ is defined by the following definition

Definition 2.1. For a given $T \in \mathcal{T}_h$ and a function $v\in V_h$, the discrete weak gradient $\nabla_d v \in RT_k(T)$ on $T$ is defined as the unique polynomial such that

where $\boldsymbol{n}$ is the unit outward normal on $\partial T$, $RT_k(T) = [Q_k(T)]^2+\bf{x}Q_k(T)$, and $\{v\}$ is defined in (6) and (7).

The weak gradient operator $\nabla_{d}$ as defined in (10) is a local operator computed at each element. It can be extended to any function $v \in V_h$ by taking weak gradient locally on each element T. More precisely, the weak gradient of any $v \in V_h$ is defined element-by-element as follows:

We introduce the following bilinear form:

the conforming DG algorithm to solve the problems (1) - (2) is given by

Conforming DG algorithm 1. Find $u_h\in V_h$ satisfying $u_h = I_hg \; \mathit{\text{on}} \; \partial \Omega$ and

where $I_h$ is the kth order Lagrange interpolation.

3.   Stability and well-posedness

We will prove the existence and uniqueness of the solution of equation (11). Firstly, we present the following two useful inequalities to derive the forthcoming analysis.

Lemma 3.1 (trace inequality). Let $T$ be an element of the finite element partition $\mathcal{T}_h$, and $e$ is an edge or face which is part of $\partial T$. For any function $\varphi\in H^1(T)$, the following trace inequality holds true (see [20] for details):

where $C$ is a constant independent of $h$.

Lemma 3.2 (inverse inequality). Let $\mathcal{T}_h$ be a finite element partition of $\Omega$ that is shape regualr. Assume that $\mathcal{T}_h$ satisfies all the assumptions A1-A4 in [20]. Then, for any piecewise polynomial function $\varphi$ of degree $n$ on $\mathcal{T}_h$, there exists a constant $C = C(n)$ such that

Then, we define the following semi-norms in the discontinuous finite element space $V_h$

We have the equivalence between the semi-norms $|\!|\!| v |\!|\!|$ and $\|v\|_{1,h}$, and it is proved in the following lemma.

Lemma 3.3. For any $v\in V_h$, the following equivalence holds true

where $C_1$ and $C_2$ are two constants independent of $h$.

Proof. It follows from the definition of $\nabla_d v$, integration by parts, the trace inequality, and the inverse inequality that

For any $e\subset\partial T_1$, $e = \partial T_1\cap\partial T_2$, we have

Then we can get

Substituting (18) into (17) gives

this completes the proof of the right-hand of (16).

To prove the left-hand of (16), we consider the subspace of $RT_k(T)$ for any $T\in\mathcal{T}_h$

Note that $D(k,T)$ is a dual space of $[Q_{k-1}(T)]^2$ [13]. Thus, for any $\nabla v\in[Q_{k-1}(T)]^2$, we have

Using the integration by parts, Cauchy-Schwarz inequality, the definition of $D(k,T)$ and $\nabla_d v$, we get

where we have used the fact that $\boldsymbol{q</italic>}\cdot \boldsymbol{n}|_{\partial T} = 0$ in the definition of $D(k,T)$. Combining the above result with (19), one has

We define the space $D_e(k,T)$ as the set of all $\boldsymbol{q</italic>} \in RT_k(T)$ such that all degrees of freedom, except those for $\boldsymbol{q</italic>} \cdot\boldsymbol{n}|_e$, vanish. Note that $D_e(k,T)$ is a dual space of $[Q_k(e)]^2$ [13]. Thus, we know

Following the integration by parts and the definition of $\nabla_d$, we can derive that

Together with (20), we obtain

Substituting the above inequality into (21), by the scaling argument [13], for such $\boldsymbol{q} \in D_e(k,T)$, we have $\|\boldsymbol{q}\|_T \leq h^{\frac{1}{2}}\|\boldsymbol{q}\cdot\boldsymbol{n}\|_e$, then

Combining (20) and (22) gives a proof of the left-hand of (16).

Lemma 3.4. The semi-norm $|\!|\!| \cdot |\!|\!|$ defined in (14) is a norm in $V_h^0$.

Proof. We shall only verify the positivity property for $|\!|\!| \cdot |\!|\!|$. To this end, assume $|\!|\!| v |\!|\!| = 0$ for some $v\in V_h^0$. By Lemma 3.3, it follows that $\| v \|_{1,h} = 0$ for all $T\in\mathcal{T}_h$, which means that $\nabla v = {\mathit{\boldsymbol{0}}}$ for all elements $T\in\mathcal{T}_h$ and $[\![ v]\! ] = 0$ for all edges $e\in\mathcal{E}_h^0$. We can derive from $\nabla v = {\mathit{\boldsymbol{0}}}$ for all $T\in\mathcal{T}_h$ that $v$ is a constant in each $T$. $[\![ v]\! ] = 0$ on each $e\in\mathcal{E}_h^0$ implies $v$ is a continuous function. This two conclusions and $v = 0$ on $\partial \Omega$ show that $v = 0$, which completes the proof of the lemma.

The above two lemmas imply the well posedness of the scheme (11). We prove the existence and uniqueness of solution of the conforming DG method in Theorem 3.1.

Theorem 3.1. The conforming DG scheme (11) has and only has one solution.

Proof. To prove the scheme (11) is uniquely solvable, it suffices to verify that the homogeneous equation has zero as its unique solution. To this end, let $u_h\in V_h$ be the solution of the numerical scheme 11 with homogeneous data $f = {0}$, ${g} = 0$. Letting $v_h = u_h$, we obtain

which leads to $u_h = 0$ by using Lemma 3.4. This completes the proof of the theorem.

4.   Error equation

In this section, we will derive an error equation which will be used for the error estimates. For any $\boldsymbol{q</italic>}\in H(div, \Omega)$, we assume that there exist an interpolation operator $\Pi_h$ satisfying $\Pi_h\boldsymbol{q}\in H(div, \Omega)\cap RT_k(T)$ on each element $T\in\mathcal{T}_h$ and

For any $w\in H^{1+k}(\Omega)$ with $k\geq 1$, from Lemma 7.3 in [20], we have the estimate of $\Pi_h$ as follows.

Moreover, it is easy to verify the following property holds true.

Lemma 4.1. For any $\boldsymbol{q}\in H(div, \Omega)$,

Proof. $\Pi_h\boldsymbol{q}\in H(div, \Omega)$ implies that $\Pi_h\boldsymbol{q}$ is continuous across each interior edge. Since $v\in V_h^0$, we know that $\{v\} = v = 0$ on $\partial \Omega$. Then

By the definition of $\Pi_h$ and $\nabla_d$ and the equation (26), we have

This completes the proof of the lemma.

Before establishing the error equation, we define a continuous finite element subspace of $V_h$ as follows

so as a subspace of $\tilde{V}_h$

Lemma 4.2. For any $v\in\tilde{V}_h$, we have

Proof. By the definition of $\nabla_d$ and integration by parts, for any $\boldsymbol{q</italic>}\in RT_k(T)$, we have

which gives

Letting $\boldsymbol{q}$ be $\nabla_d v-\nabla v$ in the above equation yields $\|\nabla_d v-\nabla v\| = 0$, which completes the proof of the lemma.

Let $e_h = I_h u - u_h$, where $I_h$ is the kth order Lagrange interpolation, $u\in H^{k+1}(\Omega)$ with $k \geq 1$ is the exact solution of the Poisson equations (1) - (2), and $u_h\in V_h$ is the numerical solution of the scheme (11). The following estimate of the Lagrange interpolation operator $I_h$ holds true.

It is obvious that $e_h\in V_h^0$ and $I_h u\in\tilde{V}_h$. We have the following lemma:

Lemma 4.3. Denote $e_h = I_h u - u_h$ the error of conforming DG method arising from (11). For any $v_h\in V_h^0$, we have

where

Proof. Since $I_h u\in\tilde{V}_h$, we have $\nabla_d I_h u = \nabla I_h u$. Using the property (25), we can derive

By the definition of the scheme (11), we have

This completes the proof of the lemma.

5.   Error estimates

The goal of this section is to derive the error estimates in $H^1$ and $L^2$ norms for the conforming DG solution $u_h$.

Theorem 5.1. Let $u\in H^{k+1}(\Omega)$ with $k \geq 1$ be the exact solution of the Poisson equation (1) - (2), and $u_h\in V_h$ be the numerical solution of the scheme (11). Let $e_h = I_h u - u_h$, there exists a constant C independent of $h$ such that

Proof. Letting $v_h = e_h$ in (31), and by the definition of $|\!|\!|\cdot|\!|\!|$, we have

From the Cauchy-Schwarz inequality, the triangle inequality, the definition of $|\!|\!|\cdot|\!|\!|$, (24), and (30), we arrive at

Then, we have

Substituting (35) to (34), we obtain

which completes the proof of the lemma.

It is obvious that $\tilde{V}_h^0\subset V_h^0$. Let $\tilde{u}_h\in\tilde{V}_h$ be the finite element solution for the problem (1)-(2) which satisfies $\tilde{u}_h = I_hg$ on $\partial\Omega$ and

For any $v\in\tilde{V}_h^0\subset\tilde{V}_h$, we have $\nabla_d v = \nabla v$, i.e.

In the rest of this section, we derive an optimal order error estimate for the conforming DG approximation (11) in $L^2$ norm by adopting the duality argument. To this end, we consider the following dual problem that seeks $\Phi\in H_0^1(\Omega)$ satisfying

Assume that the dual problem satisfies $H^2$-regularity, which means the following priori estimate holds true

In the following of this paper, we note $\varepsilon_h = u_h-\tilde{u}_h$ for simplicity.

Theorem 5.2. Assume $u\in H^{k+1}(\Omega)$ with $k \geq 1$ is the exact solution of the Poisson equation (1) - (2), and $u_h\in V_h$ is the numerical solution obtained with the scheme (11). Furthermore, assume that (39) holds true. Then, there exists a constant C independent of $h$ such that

Proof. First, we shall derive the optimal order for $\varepsilon_h$ in $L^2$ norm. Consider the corresponding conforming DG scheme defined in (11) and let $\Phi_h\in V_h^0$ be the solution satisfying

Since $I_h\Phi\in\tilde{V}_h$, it follows from (37) that

which gives

Setting $v = \varepsilon_h$ in (41), then by the definition of $\varepsilon_h$ and (42), we have

Then, by the Cauchy-Schwarz inequality, (33) and (39), we obtain

which gives

Combining the error estimate of finite element solution, the triangle inequality and (43) yields (40), which completes the proof of the theorem.

6.   Numerical experiments

In this section, we shall present some numerical results for the conforming discontinuous Galerkin method analyzed in the previous sections.

We solve the following Poisson equation on the unit square domain $\Omega = (0,1)\times (0,1)$,

The exact solution of the above problem is $u = \sin(\pi x)\sin(\pi y)$. Uniform square grids as shown in Figure 1 are used for computation.

We first use the $P_k$ conforming discontinuous Galerkin spaces (8) to compute the test case (44)-(45), where $P_k$ denotes the set of polynomials of 2 variables of degree less than or equal to k. The weak gradient is computed locally using rectangular $RT_k$ polynomials. The errors and the order of convergence of the conforming DG approximations are listed in Table 1. Optimal order of convergence is achieved in every case, which is consistent with our theory. In particular, a superconvergence of order $\mathcal{O}(h^2)$ was observed in the discrete $H^1$ norm for $P_0$ elements. Furthermore, the results obtained with $P_0$ elements seems to be slightly better than that obtained with $P_1$ elements.

The same test case is also computed using the $Q_k$ conforming DG finite element space, where $Q_k$ denotes the set of polynomials of 2 variables defined on $\Omega$, and for each variable, the degree of the variable is at most k. Table 2 illustrates the numerical performance of the corresponding conforming DG scheme. It can be seen from numerical computing that, in this case, the results obtained with the $Q_1$ element are more accurate than those obtained with $Q_0( = P_0)$ elements (see Table 1). All numerical results converge at the corresponding optimal order, which is consistent with the theory.

To test the superconvergence of $P_0$ DG element, we solve the following 2nd order elliptic equation on the unit square domain $\Omega = (0,1)\times (0,1)$,

where $f$ is chosen so that the exact solution is not symmetric,

Uniform square grids as shown in Figure 1 are used for numerical computation. The numerical results are listed in Table 3. Surprising, for this problem, the $H^1$-like norm of error superconverges at 1.5 order, and the $L^2$ error has one order of superconvergence. But we do not yet know if such a superconvergence exists in general.

To test further the superconvergence of $P_0$ DG element, we solve the following 2nd order elliptic equations on the unit square domain $\Omega = (0,1)\times (0,1)$,

where $a = 1+x+y$ and $f$ is chosen so that the exact solution is not symmetric,

Uniform square grids as shown in Figure 1 are used for computation. The numerical results are listed in Table 4. Surprising, again, the $H^1$-like norm of error superconverges at 1.5 order, and the $L^2$ error has one order of superconvergence for this problem.

7.   Conclusion

In this paper, we establish a new numerical approximation scheme based on the rectangular partition to solve second order elliptic equation. We derived the numerical scheme and then proved the optimal order of convergence of the error estimates in $L^2$ and $H^1$ norms of the conforming DG method. Numerical experiments are then present to verify the theoretical analysis, and all numerical results converging at the corresponding optimal order. Comparing with existing numerical methods, the confoming DG method has the following two characteristics: 1. The formulation is relatively simple. The stabilizer $s(\cdot\; , \; \cdot)$ is no longer needed, and the boundary function $u_b$ is omitted, which is replaced by the average of internal function $u_0$; 2. The projection operator $Q_h$ used in the traditional WG method is replaced by the Lagrange interpolation operator $I_h$, which makes the theoretical analysis much easier. As can be seen from the numerical examples in Section 6, this method reduces the programming complexity while ensuring the optimal order of convergence.