A Formulation of the Porous Medium Equation with Time-Dependent Porosity: A Priori Estimates and Regularity Results

Koffi B. Fadimba

Department of Computer Science, Engineering, and Mathematics, University of South Carolina Aiken, Aiken, USA.

DOI: 10.4236/am.2024.1510043                

Abstract

We consider a generalized form of the porous medium equation where the porosity $\varphi$ is a function of time $t$ : $\varphi =\varphi \left(x,t\right)$ : $\frac{\partial \left(\varphi S\right)}{\partial t}-\nabla \cdot \left(k\left(S\right)\nabla S\right)=Q\left(S\right).$ In many works, the porosity $\varphi$ is either assumed to be independent of (or to depend very little of) the time variable $t$ . In this work, we want to study the case where it does depend on $t$ (and $x$ as well). For this purpose, we make a change of unknown function $V=\varphi S$ in order to obtain a saturation-like (advection-diffusion) equation. A priori estimates and regularity results are established for the new equation based in part on what is known from the saturation equation, when $\varphi$ is independent of the time $t$ . These results are then extended to the full saturation equation with time-dependent porosity $\varphi =\varphi \left(x,t\right)$ . In this analysis, we make explicitly the dependence of the various constants in the estimates on the porosity $\varphi$ by the introduced transport vector $w$ , through the change of unknown function. Also we do not assume zero-flux boundary, but we carry the analysis for the case $Q\equiv 0$ .

Keywords

Porous Medium Equation, Porosity, Saturation Equation, A Priori Estimates, Regularity Results

1. Introduction

The equation:

$\frac{\partial \left(\varphi S\right)}{\partial t}+\nabla \cdot \left(f\left(S\right)u\right)-\nabla \cdot \left(k\left(S\right)\nabla S\right)=Q\left(S\right),$ (1.1)

which is often called the Saturation Equation, has been extensively studied in the literature with various simplifications (see for instance [1]-[5]). This can be seen as a generalization of the commonly studied Porous Medium Equation (PME):

$\frac{\partial S}{\partial t}-\Delta \left(S^m\right)=Q\left(S\right)$ (1.2)

(see for instance [4] [6]), where $m>1$ .

Here we want to consider first the case of (1.1) where $f\left(s\right)\equiv 0$ (or $u\equiv 0$ ), i.e. without the transport term, that is

$\frac{\partial \left(\varphi S\right)}{\partial t}-\nabla \cdot \left(k\left(S\right)\nabla S\right)=Q\left(S\right).$ (1.3)

Thus equation (1.3) is a special case of (1.1) and is a generalization of (1.2).

Equation (1.1) is obtained through modeling two-phase flow through a porous medium. Here we consider the case where there are only two phases (the wet phase and the dry phase) occupying entirely the pores of the medium, e.g. water and oil, so that the sum of the saturations of the phases equals 1. We can think of the saturation of a phase as being the “rate of presence” in the medium, in terms of volume, of the phase, as compared to the total volume of the void. The saturation equation is, in reality, coupled with a second equation relative to the Darcy law (see [7]-[12] and the bibliography therein). In this work, we consider only the saturation equation and assume that the Darcy velocity $u$ is sufficiently regular for our analysis. For other works close to the saturation problem, see, for instance, [2] [13]-[15], among others.

In the literature, equation (1.1) is considered in the form

$\varphi \frac{\partial S}{\partial t}+\nabla \cdot \left(f\left(S\right)u\right)-\nabla \cdot \left(k\left(S\right)\nabla S\right)=Q\left(S\right)$ (1.4)

either because one assumes that the porosity $\varphi$ is independent of time $t$ (or changes little with time) or because one has buried the term $S\frac{\partial \varphi }{\partial t}$ in the right hand side of (1.1). In either case, unless one keeps track of this fact, one will lose some information on the problem, especially during numerical treatment, when the porosity does vary with time.

We would like to keep the presence of the porosity $\varphi$ in the resolution of the problem. In this analysis, $\varphi$ will be present through the vector $w=\frac{\nabla \varphi }{\varphi }$ , and through ${\Vert w_t\Vert }_{\infty }$ and ${\Vert w\Vert }_{\infty }$ . So, for the analysis. we assume that all second order partial derivatives of $\varphi$ exist and are bounded on $\text{Ω}\times \left(0,T_0\right]$ .

If, in (1.4), we make the additional assumption that $u=0$ and $\varphi =1$ , then we obtain

$\frac{\partial S}{\partial t}-\nabla \cdot \left(k\left(S\right)\nabla S\right)=Q\left(S\right),$ (1.5)

so that the classical porous medium equation (1.2) is a special case of (1.5) with

$k\left(s\right)=m{s}^{m-1}$

in that case.

In this paper, we investigate the case where $\varphi$ is allowed to change with time and space:

$\varphi =\varphi \left(x,t\right).$ (1.6)

For an example of time-dependent porosity, we refer to [16] where the author considers the case where the porosity $\varphi$ is a function of time $t$ (but not of the spatial variable $x$ ) for a one-dimensional Boussinesq equation for flow in an aquifer. See also [17].

We will first study a case of (1.1) where there is no transport term i.e. we want to study the generalized porous medium problem

$\frac{\partial \left(\varphi S\right)}{\partial t}-\nabla \cdot \left(k\left(S\right)\nabla S\right)=Q\left(S\right)\text{in}\Omega \times \left(0,T_0\right)$ (1.7)

$-k\left(S\right)\frac{\partial S}{\partial n}=g_1\text{on}\partial \Omega \times \left(0,T_0\right)$ (1.8)

$S\left(x,0\right)=S_0\left(x\right)\text{in}\Omega ,$ (1.9)

where Ω is a bounded domain of $R^n$ , $n=1,2,3$ , and $T_0>0$ .

Unlike the case of the classical porous medium equation, we assume that

$k\left(0\right)=k\left(1\right)=0,$ (1.10)

to account for the fact that, in this work, we have in mind an immiscible two-phase flow (see, for instance, [7] [9] [10] [18]). We assume, for the purpose of this analysis, that the porosity $\varphi$ of the medium satisfies the condition

$0<{\varphi }_0\le \varphi \left(x,t\right)\le {\varphi }_1<\infty \text{for}\text{all}\left(x,t\right)\in \Omega \times \left[0,T_0\right].$ (1.11)

For the needs of the analysis in this paper, unless otherwise explicitly stated, we make the following assumptions on $k$ .

$k\left(\xi \right)\ge \{\begin{array}{ll}{c}_1{\xi }^{\mu }\hfill & \text{if}0\le \xi \le {\alpha }_1\hfill \\ c_2\hfill & \text{if}{\alpha }_1<\xi <{\alpha }_2\hfill \\ c_3{\left(1-\xi \right)}^{\mu }\hfill & \text{if}{\alpha }_2\le \xi \le 1\hfill \end{array}$ (1.12)

where $0<{\alpha }_1<\frac{1}{2}<{\alpha }_2<1$ , and $0<\mu \le 2$ .

This work will follow, in part, the layout of a chapter contribution by the author in [19] with improved proofs and additional results.

This contribution is structured as follows. In section 2, the problem is reformulated through a change of the unknown function. For completeness, we state and establish some preliminary facts about the data. In particular, we revisit a lemma which was proved in [2], for $\varphi =1$ , and reprove it for any $\varphi$ , for the use in the context of this work. In section 3, we establish improved a priori estimates for the new problem, which include non-homogeneous boundary conditions as well as some regularity results for the solution of the new problem, with Lemma 3.7 being new at our knowledge. In section 4, we go back to the full saturation problem (1.1) and establish a priori estimates and regularity results for the case when the porosity $\varphi$ is not necessarily independent of the time variable $t$ .

Standard mathematical notations are used in this paper. For example, we define $\left(f,g\right):={\left(f,g\right)}_{\Omega }:={\displaystyle {\int }_{\Omega }fg\text{d}x}$ when this has a meaning, and ${\Vert f\Vert }_{\infty }:={\Vert f\Vert }_{L^{\infty }\left(0,T_0;L^{\infty }\left(\text{Ω}\right)\right)}$ .

2. Reformulation of the Problem

We make a change of unknown function that absorbs the porosity $\varphi$ in the first term of (1.4) i.e. we make the transformation

$V=\varphi S,$ (2.1)

(see also [16] for the case of one dimensional flow in an aquifer). Because of the assumption (1.11), we can always get back the saturation $S$ by dividing by $\varphi$ . Our approach has at least the advantage of making explicit the dependence of constants, in the analysis, on the porosity $\varphi$ , both as a function of the temporal

variable $t$ and the spatial variable $x$ , through $w=\frac{\nabla \varphi }{\varphi }$ and $\frac{\partial w}{\partial t}$ , for instance.

Using this change of unknown function, (1.7) becomes

$\frac{\partial V}{\partial t}-\nabla \cdot \left(k\left(\frac{V}{\varphi }\right)\nabla \left(\frac{V}{\varphi }\right)\right)=Q\left(\frac{V}{\varphi }\right).$ (2.2)

We define the new functions:

$\tilde{k}\left(v\left(\cdot ,\cdot \right)\right):=\frac{1}{\varphi \left(\cdot ,\cdot \right)}k\left(\frac{v\left(\cdot ,\cdot \right)}{\varphi \left(\cdot ,\cdot \right)}\right)=\frac{1}{\varphi \left(\cdot ,\cdot \right)}k\left(s\left(\cdot ,\cdot \right)\right),$ (2.3)

with $v=\varphi s$ , and

$\tilde{Q}\left(v\left(\cdot ,\cdot \right)\right):=Q\left(\frac{v\left(\cdot ,\cdot \right)}{\varphi \left(\cdot ,\cdot \right)}\right)=Q\left(s\left(\cdot ,\cdot \right)\right).$ (2.4)

Now notice that

$\nabla \left(\frac{V}{\varphi }\right)=\frac{\nabla V\varphi -V\nabla \varphi }{{\varphi }^2}.$ (2.5)

Substituting in (2.2) and using (2.4), we obtain

$\frac{\partial V}{\partial t}-\nabla \cdot \left(\tilde{k}\left(V\right)\nabla V-\tilde{k}\left(V\right)\frac{V\nabla \varphi }{\varphi }\right)=\tilde{Q}\left(V\right).$ (2.6)

With this change of unknown function we can rewrite (1.7) as follows,

$\frac{\partial V}{\partial t}+\nabla \cdot \left(\tilde{k}\left(V\right)V\frac{\nabla \varphi }{\varphi }\right)-\nabla \cdot \left(\tilde{k}\left(V\right)\nabla V\right)=\tilde{Q}\left(V\right).$ (2.7)

Next, let

$D\left(v\right)=\tilde{k}\left(v\right)$ (2.8)

$F\left(v\right)=\tilde{k}\left(v\right)v=D\left(v\right)v$ (2.9)

$w=\frac{\nabla \varphi }{\varphi },$ (2.10)

then (2.7) takes the advection-diffusion equation form:

$\frac{\partial V}{\partial t}+\nabla \cdot \left(F\left(V\right)w\right)-\nabla \cdot \left(D\left(V\right)\nabla V\right)=\tilde{Q}\left(V\right),$ (2.11)

except that the function $F$ does not satisfy exactly the same conditions as the fractional flow function from the saturation equation (1.1). Next, we have

$\begin{array}{c}\left(D\left(V\right)\nabla V\right)\cdot n=\frac{1}{\varphi }k\left(S\right)\nabla \left(\varphi S\right)\cdot n\\ =\frac{1}{\varphi }k\left(S\right)\left(S\nabla \varphi +\varphi \nabla S\right)\cdot n\\ =-g_1+F\left(V\right)w\cdot n,\end{array}$ (2.12)

when we use (1.8), (2.3), (2.9), and (2.10). Hence, we obtain the new advection-diffusion problem:

$\{\begin{array}{ll}\frac{\partial V}{\partial t}+\nabla \cdot \left(F\left(V\right)w\right)-\nabla \cdot \left(D\left(V\right)\nabla V\right)=\tilde{Q}\left(V\right)\hfill & \text{on}\Omega \times \left(0,T_0\right]\hfill \\ \left(F\left(V\right)w-\left(D\left(V\right)\nabla V\right)\right)\cdot n=g_1\hfill & \text{on}\partial \Omega \times \left[0,T_0\right]\hfill \\ V\left(x,0\right)=V_0\left(x\right)=\varphi \left(x,0\right)S_0\left(x\right)\hfill & \text{on}\Omega \hfill \end{array}$ (2.13)

We assume that:

$g_1\in L^{\infty }\left(0,T_0,L^{\infty }\left(\partial \Omega \right)\right).$ (2.14)

Remark 2.1 Even though the porosity $\varphi$ is buried in the variable $V$ , its presence in the equation is explicitly signaled by the vector $w=\frac{\nabla \varphi }{\varphi }$ .

One can easily check that $D$ verifies the following conditions.

$D\left(0\right)=D\left(\varphi \right)=0,$ (2.15)

$D\left(\xi \right)\ge \{\begin{array}{ll}{c}_1{\xi }^{\mu }\hfill & \text{if}0\le \xi \le {\alpha }_1\varphi \hfill \\ c_2\hfill & \text{if}{\alpha }_1\varphi <\xi <{\alpha }_2\varphi \hfill \\ c_3{\left(\varphi -\xi \right)}^{\mu }\hfill & \text{if}{\alpha }_2\varphi \le \xi \le \varphi \hfill \end{array}$ (2.16)

where $0<{\alpha }_1<\frac{1}{2}<{\alpha }_2<1$ , and $0<\mu \le 2$ .

For example,

$D\left(\xi \right)=\frac{1}{\varphi }k\left(\frac{\xi }{\varphi }\right)\ge c_3\frac{1}{\varphi }{\left(1-\frac{\xi }{\varphi }\right)}^{\mu }=c_3\frac{1}{{\varphi }^{\mu +1}}{\left(\varphi -\xi \right)}^{\mu }$ (2.17)

if ${\alpha }_2\le \frac{\xi }{\varphi }\le 1$ , i.e. if ${\alpha }_2\varphi \le \xi \le \varphi .$ Now, by (1.11), $\frac{c_3}{{\varphi }^{\mu +1}}\ge C_3$ , for some positive constant $C_3$ , which we still call $c_3$ . The other inequalities in (2.16) are obtained in the same manner.

We clearly have:

$F^{\prime }\left(v\right)=\frac{1}{{\varphi }^2}{k}^{\prime }\left(\frac{v}{\varphi }\right)v+\frac{1}{\varphi }k\left(\frac{v}{\varphi }\right).$ (2.18)

So, we see that $F^{\prime }\left(0\right)=F^{\prime }\left(\varphi \right)=0$ , when we consider (1.10) and (4.10).

Finally, define $\mathcal{K}$ by

$\mathcal{K}\left(v\right)={\displaystyle {\int }_0^{v}D\left(\tau \right)\text{d}\tau }.$ (2.19)

For the remainder of this paper, in order to simplify the analysis, we assume the following, unless otherwise explicitly stated.

$\left|\Omega \right|=1.$ (2.20)

Therefore, the focus of this work is the problem:

$\{\begin{array}{ll}\frac{\partial V}{\partial t}+\nabla \cdot \left(F\left(V\right)w\right)-\Delta \mathcal{K}\left(V\right)=\tilde{Q}\left(V\right)\hfill & \text{on}\Omega \times \left(0,T_0\right]\hfill \\ \left(F\left(V\right)w-\nabla \left(\mathcal{K}\left(S\right)\right)\right)\cdot n=g_1\hfill & \text{on}\partial \Omega \times \left[0,T_0\right]\hfill \\ V\left(x,0\right)=V_0\left(x\right):=\varphi \left(x,0\right)S_0\left(x\right)\hfill & \text{on}\Omega \hfill \end{array}$ (2.21)

Under conditions (2.16), (1.10), and (4.10), we have the following lemma.

Lemma 2.2 If the function $F$ is twice continuously differentiable near 0 and $\varphi$ , then

$|F\left(v_2\right)-F\left(v_1\right){|}^2\le C^{\text{*}}\left(\mathcal{K}\left(v_2\right)-\mathcal{K}\left(v_1\right)\right)\left(v_2-v_1\right)$ (2.22)

for all $0\le v_1\le v_2\le \varphi$ , where $C^{\text{*}}=C^{\text{*}}\left(\mu ,\varphi \right)>0$ .

Lemma 2.2 has been proved in [2] for $\varphi \equiv 1$ (the assumption $\varphi \equiv 1$ is often made for mathematical analysis purpose only), but here we give a proof for any $\varphi$ and for ${\alpha }_2\varphi \le v_1\le v_2\le \varphi$ only. For $0\le v_1\le v_2\le {\alpha }_1\varphi \le 1$ , the proof goes the same way: see [19]. The constants ${\alpha }_1$ and ${\alpha }_2$ are as in (2.16).

Because $\mathcal{K}$ is obviously Lipschtitz (and monotone increasing, since $D\left(v\right)\ge 0$ , for all $0\le v\le \varphi$ ), the following holds:

${\left|\mathcal{K}\left(v_2\right)-\mathcal{K}\left(v_1\right)\right|}^2\le C^{\text{**}}\left(\mathcal{K}\left(v_2\right)-\mathcal{K}\left(v_1\right)\right)\left(v_2-v_1\right),$ (2.23)

for some $C^{\text{**}}>0$ and for all $0\le v_1,v_2\le \varphi$ .

Remark 2.3 Lemma 2.2 implies

$\left|F^{\prime }\left(v\right)\right|\le C^{\text{*}}\sqrt{D\left(v\right)}.$ (2.24)

for all $0\le v\le \varphi$ .

We also state the following result which was proved in in [2], based solely on condition (2.16), for the case $\varphi \equiv 1$ .

Lemma 2.4 Under condition (2.16), we have

${\left(v_2-v_1\right)}^{\mu +2}\le C\left(\mathcal{K}\left(v_2\right)-\mathcal{K}\left(v_1\right)\right)\left(v_2-v_1\right),$ (2.25)

for all $0\le v_1,v_2\le \varphi$ .

3. A Priori Estimates and Regularities Results

For a priori estimates (uniqueness) and some regularity results for Problem 2.21, we need the Poisson Solution Operator $T$ and some of its properties. We state the properties we will need in this section and other sections to follow. For proofs of these properties and more details see, for instance, [1] [2] [20].

3.1. The Poisson Solution Operator

Consider the elliptic boundary value problem:

$\{\begin{array}{ll}-\Delta \omega =f-f_{\Omega }\hfill & \text{in}\Omega \hfill \\ \frac{\partial \omega }{\partial n}=0\hfill & \text{on}\partial \Omega \hfill \\ {\omega }_{\Omega }=f_{\Omega }\hfill & \hfill \end{array}$ (3.1)

with $f\in {\left(H^1\left(\Omega \right)\right)}^{*}$ . Then (see [21] [22]) problem 3.1 has a unique solution $\omega \in H^1$ .

We define the solution operator $T:{\left(H^1\right)}^{*}\to H^1$ by $T\left(f\right)=\omega$ , where $\omega \in H^1$ is the unique solution to (3.1), and $f\in {\left(H^1\right)}^{\text{*}}$ .

A weak formulation of (3.1) is expressed as:

$\left(\nabla \left(Tf\right),\nabla \varphi \right)=\left(f,\varphi \right)-f_{\Omega }{\varphi }_{\Omega },$ (3.2)

for all $\varphi \in H^1\left(\Omega \right)$ .

This operator satisfies the following.

Proposition 3.1 The operator $T$ , defined by (3.1), is linear, symmetric and positive definite.

$\left(Tf,g\right)=\left(f,Tg\right)\text{for}\text{all}f,g\in {\left(H^1\right)}^{\text{*}}.$ (3.3)

We also have

$\left(Tf,f\right)={\Vert \nabla Tf\Vert }_{L^2}^2+{\left(Tf\right)}_{\Omega }^2={\Vert \nabla Tf\Vert }_{L^2}^2+{\left(f\right)}_{\Omega }^2$ (3.4)

With this proposition in mind, we can define on ${\left(H^1\right)}^{\text{*}}$ the norm:

${\Vert \left|f\right|\Vert }_{{\left(H^1\right)}^{\text{*}}}:={\left(Tf,f\right)}^{\frac{1}{2}}$ (3.5)

This is the same as the dual norm for $H^1$ , when $H^1$ is equipped with the norm

${\Vert \left|\omega \right|\Vert }_{H^1}:={\left({\Vert \nabla \omega \Vert }_{L^2}^2+{\left({\omega }_{\Omega }\right)}^2\right)}^{\frac{1}{2}}.$ (3.6)

With these definitions we get the identity

$\left(Tf,f\right)={\Vert \left|Tf\right|\Vert }_{H^1}^2\simeq {\Vert f\Vert }_{{\left(H^1\right)}^{\text{*}}}^2.$ (3.7)

Proposition 3.2 Suppose f belongs to ${\left(H^1\right)}^{*}$ , then

${\left(Tf,f\right)}^{\frac{1}{2}}={\Vert f\Vert }_{{\left(H^1\right)}^{\text{*}}},$ (3.8)

in the sense of the norm (3.6).

3.2. A Priori Estimates

A weak formulation of (2.21) is:

$\left(\frac{\partial V}{\partial t},\chi \right)-\left(F\left(V\right)w,\nabla \chi \right)+\left(\nabla \mathcal{D}\left(V\right),\nabla \chi \right)-{\displaystyle {\int }_{\partial \Omega }{g}_1\chi \text{d}\sigma }=\left(\tilde{Q}\left(V\right),\chi \right),$ (3.9)

for all $\chi \in H^1\left(\Omega \right)$ and $t\in \left(0,T_0\right]$ .

We prove the theorem for any $\varphi$ , not just for $\varphi \equiv 1$ . Hence the theorem is more general than its counterpart in [19]. To simplify the analysis, we will assume that $Q\equiv 0$ , unless otherwise stated. For a more general case, with $Q\in L^{\infty }\left(0,T_0,L^2\left(\Omega \right)\right)$ , we can proceed as in [2].

Theorem 3.3 Let $V_1$ and $V_2$ be (weak) solutions to Problem 2.21, corresponding to the initial conditions $V_1^0$ and $V_2^0$ , respectively. Let $g_{1,1}$ and $g_{1,2}$ be the respective righthand sides of the corresponding boundary conditions. Then

$\begin{array}{l}\underset{0\le t\le T_0}{\sup }{\Vert V_2-V_1\Vert }_{{\left(H^1\right)}^{*}}^2+{\displaystyle {\int }_0^{T_0}\left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right),V_2-V_1\right)\left(\tau \right)\text{d}\tau }\\ \le C\max \left(1,{\Vert w\Vert }_{L^{\infty }\left(L^{\infty }\right)}\right)\left\{{\Vert V_2^0-V_1^0\Vert }_{{(H^1)}^{*}}^2+{\Vert g_{1,2}-g_{1,1}\Vert }_{L^2\left(L^2\left(\partial \Omega \right)\right)}^2\right\}.\end{array}$ (3.10)

Proof. We subtract the corresponding equation (3.9) for the triplet $\left(V_1,V_1^0,g_{1,1}\right)$ , from the equation corresponding to the triplet $\left(V_2,V_2^0,g_{1,2}\right)$ , and let $\chi =T\left(V_2-V_1\right)$ to obtain

$\begin{array}{l}\left(\frac{\partial \left(V_2-V_1\right)}{\partial t},T\left(V_2-V_1\right)\right)-\left(\left(F\left(V_2\right)-F\left(V_1\right)\right)w,\nabla T\left(V_2-V_1\right)\right)\\ +\left(\nabla \left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\right),\nabla T\left(V_2-V_1\right)\right)\\ ={\displaystyle {\int }_{\partial \Omega }\left(g_{1,2}-g_{1,1}\right)T\left(V_2-V_1\right)\text{d}\sigma }.\end{array}$ (3.11)

Using the symmetry of the operator $T$ and (3.8), we have

$\begin{array}{c}\left(\frac{\partial \left(V_2-V_1\right)}{\partial t},T\left(V_2-V_1\right)\right)=\left(T\left(\frac{\partial \left(V_2-V_1\right)}{\partial t}\right),V_2-V_1\right)\\ =\left(\frac{\partial T\left(V_2-V_1\right)}{\partial t},V_2-V_1\right)\\ =\frac{1}{2}\frac{\text{d}}{\text{d}t}\left(T\left(V_2-V_1\right),V_2-V_1\right)\\ =\frac{1}{2}\frac{\text{d}}{\text{d}t}{\Vert V_2-V_1\Vert }_{{\left(H^1\right)}^{\text{*}}}^2.\end{array}$ (3.12)

Using (3.12), we can rewrite (3.11) as follows.

$\begin{array}{l}\frac{1}{2}\frac{\text{d}}{\text{d}t}{\Vert V_2-V_1\Vert }_{{\left(H^1\right)}^{*}}^2+\left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right),V_2-V_1\right)\\ =\left(\left(F\left(V_2\right)-F\left(V_1\right)\right)w,\nabla T\left(V_2-V_1\right)\right)+{\left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\right)}_{\Omega }{\left(V_2-V_1\right)}_{\Omega }\\ +{\displaystyle {\int }_{\partial \Omega }\left(g_{1,2}-g_{1,1}\right)T\left(V_2-V_1\right) \text{d}\sigma },\end{array}$ (3.13)

where we have made use of (3.2). Now, we treat each of the two terms on the right hand side of (3.13) separately. For the first term we obtain:

$\begin{array}{l}\left|\left(\left(F\left(V_2\right)-F\left(V_1\right)\right)w,\nabla T\left(V_2-V_1\right)\right)\right|\\ \le \frac{1}{4C^{*}}{\Vert F\left(V_2\right)-F\left(V_1\right)\Vert }_{L^2}^2+C{\Vert w\Vert }_{L^{\infty }\left(L^{\infty }\right)}^2{\Vert \nabla T\left(V_2-V_1\right)\Vert }_{L^2}^2\\ \le \frac{1}{4C^{*}}{\Vert F\left(V_2\right)-F\left(V_1\right)\Vert }_{L^2}^2+C{\Vert w\Vert }_{L^{\infty }\left(L^{\infty }\right)}^2{\Vert V_2-V_1\Vert }_{{\left(H^1\right)}^{*}}^2,\end{array}$ (3.14)

by the arithmetic-geometric mean inequality. We have also used the fact that (3.4) and (3.8) imply

${\Vert \nabla T\left(V_2-V_1\right)\Vert }_{L^2}\le {\Vert V_2-V_1\Vert }_{{\left(H^1\right)}^{\text{*}}}.$ (3.15)

As for the second term, we have

$\begin{array}{c}\left|{\left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\right)}_{\Omega }\right|=\left|{\displaystyle {\int }_{\Omega }\left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\right)\text{d}x}\right|\\ \le {\Vert \mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\Vert }_{L^2}{\Vert 1\Vert }_{L^2}\\ \le {\Vert \mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\Vert }_{L^2},\end{array}$ (3.16)

by (2.20).

In the same way, one obtains:

$\left|{\left(V_2-V_1\right)}_{\Omega }\right|=\left|{\displaystyle {\int }_{\Omega }\left(V_2-V_1\right)\text{d}x}\right|\le {\Vert V_2-V_1\Vert }_{{\left(H^1\right)}^{*}}{\Vert 1\Vert }_{H^1}\le C{\Vert V_2-V_1\Vert }_{{\left(H^1\right)}^{*}},$ (3.17)

with the help of the Holder's Inequality.

Therefore, with the help of the arithmetic-geometric mean inequality, the second term on the right hand side of (3.13) can be bounded as follows.

$\begin{array}{l}\left|{\left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\right)}_{\Omega }{\left(V_2-V_1\right)}_{\Omega }\right|\\ \le \frac{1}{4C^{**}}{\Vert \mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\Vert }_{L^2}^2+C{\Vert V_2-V_1\Vert }_{{\left(H^1\right)}^{*}}^2.\end{array}$ (3.18)

We treat the last term of (3.13) in the following manner.

$\begin{array}{c}\left|{\displaystyle {\int }_{\partial \Omega }\left(g_{1,2}-g_{1,1}\right)T\left(V_2-V_1\right)\text{d}\sigma }\right|\le {\Vert g_{1,2}-g_{1,1}\Vert }_{L^2\left(\partial \Omega \right)}{\Vert T\left(V_2-V_1\right)\Vert }_{L^2\left(\partial \Omega \right)}\\ \le \frac{1}{2}{\Vert g_{1,2}-g_{1,1}\Vert }_{L^2\left(\partial \Omega \right)}^2+\frac{1}{2}{\Vert T\left(V_2-V_1\right)\Vert }_{L^2\left(\partial \Omega \right)}^2\\ \le \frac{1}{2}{\Vert g_{1,2}-g_{1,1}\Vert }_{L^2\left(\partial \Omega \right)}^2+C{\Vert T\left(V_2-V_1\right)\Vert }_{H^1\left(\Omega \right)}^2\\ \le \frac{1}{2}{\Vert g_{1,2}-g_{1,1}\Vert }_{L^2\left(\partial \Omega \right)}^2+C{\Vert V_2-V_1\Vert }_{{\left(H^1\right)}^{*}}^2,\end{array}$ (3.19)

where we have used the Trace Theorem (see [23]) followed by (3.7).

Putting (3.13), (3.14), and (3.18) together, we obtain

$\begin{array}{l}\frac{1}{2}\frac{\text{d}}{\text{d}t}{\Vert V_2-V_1\Vert }_{{\left(H^1\right)}^{*}}^2+\left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right),V_2-V_1\right)\\ \le \frac{1}{4C^{*}}{\Vert F\left(V_2\right)-F\left(V_1\right)\Vert }_{L^2}^2+\frac{1}{4C^{**}}{\Vert \mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\Vert }_{L^2}^2\\ +C{\Vert w\Vert }_{L^{\infty }\left(L^{\infty }\right)}^2{\Vert V_2-V_1\Vert }_{{\left(H^1\right)}^{*}}^2+\frac{1}{2}{\Vert g_{1,2}-g_{1,1}\Vert }_{L^2\left(\partial \Omega \right)}^2.\end{array}$ (3.20)

Finally, we hide the first two terms on the right hand side of (3.20) into the left hand sides, thanks to (2.22) and (2.23), and then apply the Gronwall Lemma to end the proof of the theorem.

The following corollary is an immediate consequence of Theorem 3.3.

Corollary 3.4 Under the conditions of Theorem 3.3, we have

${\Vert \mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\Vert }_{L^2\left(L^2\right)}\le C\left\{{\Vert V_2^0-V_1^0\Vert }_{{\left(H^1\right)}^{\text{*}}}+{\Vert g_{1,2}-g_{1,1}\Vert }_{L^2\left(L^2\left(\partial \Omega \right)\right)}\right\},$ (3.21)

${\Vert V_2-V_1\Vert }_{L^{2+\mu }\left(L^{2+\mu }\right)}^{2+\mu }\le C\left\{{\Vert V_2^0-V_1^0\Vert }_{{\left(H^1\right)}^{*}}^2+{\Vert g_{1,2}-g_{1,1}\Vert }_{L^2\left(L^2\left(\partial \Omega \right)\right)}^2\right\},$ (3.22)

and

${\Vert V_2-V_1\Vert }_{L^2\left(L^2\right)}\le C\left\{{\Vert V_2^0-V_1^0\Vert }_{{\left(H^1\right)}^{*}}^{\frac{2}{2+\mu }}+{\Vert g_{1,2}-g_{1,1}\Vert }_{L^2\left(L^2\left(\partial \Omega \right)\right)}^{\frac{2}{\mu +2}}\right\}.$ (3.23)

The following two regularity results help establish a priori estimates that were not established in [2] for the case $\varphi \equiv 1$ .

The following lemma gives regularities on the solution $V$ of (2.21).

Lemma 3.5 Let $V$ be the solution to Problem 2.21. Then

${\Vert V\Vert }_{L^{\infty }\left(L^2\right)}^2+\sqrt{D\left(V\right)}\nabla V_{L^2\left(L^2\right)}^2\le C,$ (3.24)

where $C=C\left({\Vert w\Vert }_{\infty },T_0,{\Vert V^0\Vert }_{L^2},{\Vert g_1\Vert }_{\infty }\right)$ .

Proof. If we let $\chi =V$ in (3.9) we obtain, in a distributional sense,

$\left(\frac{\partial V}{\partial t},V\right)-\left(F\left(V\right)w,\nabla V\right)+\left(\nabla \mathcal{K}\left(V\right),\nabla V\right)-{\displaystyle {\int }_{\partial \Omega }{g}_{1}V\text{d}\sigma }=0,$ (3.25)

which can be rewritten as follows.

$\frac{1}{2}\frac{\text{d}}{\text{d}t}{\Vert V\Vert }_{L^2}^2+\sqrt{D\left(V\right)}\nabla V_{L^2}^2=\left(F\left(V\right)w,\nabla V\right)+{\displaystyle {\int }_{\partial \Omega }{g}_{1}V\text{d}\sigma },$ (3.26)

since ${\mathcal{K}}^{\prime }\left(V\right)=D\left(V\right)$ by (2.19). Now $F\left(V\right)=D\left(V\right)V$ , by (2.9) and (2.8), so that we get, by (3.26),

$\begin{array}{l}\frac{1}{2}\frac{\text{d}}{\text{d}t}{\Vert V\Vert }_{L^2}^2+{\Vert \sqrt{D\left(V\right)}\nabla V\Vert }_{L^2}^2=\left(Vw,D\left(V\right)\nabla V\right)+{\displaystyle {\int }_{\partial \Omega }{g}_{1}V\text{d}\sigma }\\ \le \frac{1}{2}{\Vert \sqrt{D\left(V\right)}Vw\Vert }_{L^2}^2+\frac{1}{2}{\Vert \sqrt{D\left(V\right)}\nabla V\Vert }_{L^2}^2+{\Vert g_{1}V\Vert }_{\infty }\left({\text{vol}}_{n-1}\left(\partial \Omega \right)\right)\\ \le {\Vert \sqrt{D\left(V\right)}w\Vert }_{L^{\infty }\left(L^{\infty }\right)}^2{\Vert V\Vert }_{L^2}^2+\frac{1}{2}{\Vert \sqrt{D\left(V\right)}\nabla V\Vert }_{L^2}^2+{\Vert g_{1}V\Vert }_{\infty }\left({\text{vol}}_{n-1}\left(\partial \Omega \right)\right)\\ \le C\left({\Vert w\Vert }_{L^{\infty }\left(L^{\infty }\right)}\right){\Vert V\Vert }_{L^2}^2+\frac{1}{2}{\Vert \sqrt{D\left(V\right)}\nabla V\Vert }_{L^2}^2+{\Vert g_{1}V\Vert }_{\infty }\left({\text{vol}}_{n-1}\left(\partial \Omega \right)\right).\end{array}$ (3.27)

We hide the second term on the right hand side of (3.27) into the like term on the left hand side, Apply the Gronwall Lemma over the interval $\left[0,t\right]$ , $0<t\le T_0$ , and then take the supremum over $\left[0,T_0\right]$ , to get

$\begin{array}{l}{\Vert V\Vert }_{L^{\infty }\left(L^2\right)}^2+{\Vert \sqrt{D\left(V\right)}\nabla V\Vert }_{L^2\left(L^2\right)}^2\\ \le C\left({\Vert w\Vert }_{L^{\infty }\left(L^{\infty }\right)},T_0,{\Vert V^0\Vert }_{L^2}^2\right)+C{\Vert g_1\Vert }_{\infty }\left({\text{vol}}_{n-1}\left(\partial \Omega \right)\right).\end{array}$ (3.28)

Hence (3.24).

Remark 3.6

1. Lemma 3.5 states that, if $V$ is a weak solution of Problem 2.21, then

$V\in L^{\infty }\left(0,T_0;L^2\left(\Omega \right)\right),$ (3.29)

and

$\mathcal{K}\left(V\right)\in L^2\left(0,T_0;H^1\left(\Omega \right)\right).$ (3.30)

The latter is true because

$\begin{array}{c}{\Vert \nabla \mathcal{K}\left(V\right)\Vert }_{L^2}={\Vert D\left(V\right)\nabla V\Vert }_{L^2}={\Vert \sqrt{D(\cdot )\left(V\right)}\Vert }_{\infty }{\Vert \sqrt{D\left(V\right)}\nabla V\Vert }_{L^2}\\ \le C{\Vert \sqrt{D\left(V\right)}\nabla V\Vert }_{L^2}.\end{array}$ (3.31)

2. Notice that, in the proof of Lemma 3.5, condition (2.16) was not used. This point will be of great importance later when we tackle the proof of existence and uniqueness for Problem 2.21, in a sequel of this paper.

The next lemma proves the following regularity results for $\left(\mathcal{K}\left(V\right)\right)$ and ${\left(\mathcal{K}\left(V\right)\right)}_t$ .

$\mathcal{K}\left(V\right)\in L^{\infty }\left(0,T_0;H^1\left(\Omega \right)\right)$ (3.32)

and

${\left(\mathcal{K}\left(V\right)\right)}_t\in L^2\left(0,T_0;L^2\left(\Omega \right)\right).$ (3.33)

Lemma 3.7 Let $V$ be a weak solution of Problem 2.21. Then

${\Vert \sqrt{D\left(V\right)}{V}_t\Vert }_{L^2\left(L^2\right)}^2+\eta {\Vert \nabla \mathcal{K}\left(V\right)\Vert }_{L^{\infty }\left(L^2\right)}^2\le C.$ (3.34)

where $\eta >0$ and $C=C\left(V^0,{\Vert w\Vert }_{L^{\infty }\left(L^{\infty }\right)},{\Vert w_t\Vert }_{L^2\left(L^2\right)},{\Vert g_1\Vert }_{\infty },{\Vert {\left(g_1\right)}_t\Vert }_{\infty }\right)$ .

Proof. Set $\chi ={\left(\mathcal{K}\left(V\right)\right)}_t$ in (3.9) to get

$\begin{array}{l}\left(V_t,{\left(\mathcal{K}\left(V\right)\right)}_t\right)-\left(F\left(V\right)w,\nabla {\left(\mathcal{K}\left(V\right)\right)}_t\right)+\left(\nabla \mathcal{K}\left(V\right),\nabla {\left(\mathcal{K}\left(V\right)\right)}_t\right)\\ -{\displaystyle {\int }_{\partial \Omega }{g}_1{\left(\mathcal{K}\left(V\right)\right)}_t\text{d}\sigma }=0.\end{array}$ (3.35)

Rewriting (3.35), we get

$\begin{array}{l}{\Vert \sqrt{D\left(V\right)}{V}_t\Vert }_{L^2}^2+\frac{1}{2}\frac{\text{d}}{\text{d}t}{\Vert \nabla \mathcal{K}\left(V\right)\Vert }_{L^2}^2\\ ={\left(F\left(V\right)w\cdot \nabla \left(\mathcal{K}\left(V\right)\right)\right)}_t+{\displaystyle {\int }_{\partial \Omega }{g}_1{\left(\mathcal{K}\left(V\right)\right)}_t\text{d}\sigma }.\end{array}$ (3.36)

By the product rule,

$\left(F\left(V\right)w,\nabla {\left(\mathcal{K}\left(V\right)\right)}_t\right)=\frac{\text{d}}{\text{d}t}\left(F\left(V\right)w,\nabla \left(\mathcal{K}\left(V\right)\right)\right)-\left({\left(F\left(V\right)w\right)}_t,\nabla \left(\mathcal{K}\left(V\right)\right)\right).$ (3.37)

The last term on the right hand side of (3.37) is rewritten as

$\left({\left(F\left(V\right)w\right)}_t,\nabla \left(\mathcal{K}\left(V\right)\right)\right)=\left(F^{\prime }\left(V\right)V_{t}w+F\left(V\right)w_t,\nabla \left(\mathcal{K}\left(V\right)\right)\right).$ (3.38)

Therefore (3.37) becomes

$\begin{array}{l}{\left(F\left(V\right)w,\nabla \left(\mathcal{K}\left(V\right)\right)\right)}_t\\ \le \frac{\text{d}}{\text{d}t}\left(F\left(V\right)w,\nabla \left(\mathcal{K}\left(V\right)\right)\right)+\left|\left(F^{\prime }\left(V\right)V_{t}w,\nabla \left(\mathcal{K}\left(V\right)\right)\right)\right|\\ +\left|\left(F\left(V\right)w_t,\nabla \left(\mathcal{K}\left(V\right)\right)\right)\right|\\ \le \frac{\text{d}}{\text{d}t}\left(F\left(V\right)w,\nabla \left(\mathcal{K}\left(V\right)\right)\right)+\frac{1}{2C^{*}}{\Vert \frac{F^{\prime }\left(V\right)}{\sqrt{D\left(V\right)}}\Vert }_{L^{\infty }\left(L^{\infty }\right)}^2{\Vert \sqrt{D\left(V\right)}{V}_t\Vert }_{L^2}^2\\ +C{\Vert F\left(V\right)w_t\Vert }_{L^2}^2+C\left({\Vert w\Vert }_{L^{\infty }\left(L^{\infty }\right)}\right){\Vert \nabla \left(\mathcal{K}\left(V\right)\right)\Vert }_{L^2}^2,\end{array}$ (3.39)

where $C^{\text{*}}$ is as in (2.24).

The boundary term in (3.36) is handled as follows.

${\displaystyle {\int }_{\partial \Omega }{g}_1{\left(\mathcal{K}\left(V\right)\right)}_t\text{d}\sigma }=\frac{\text{d}}{\text{d}t}\left({\displaystyle {\int }_{\partial \Omega }{g}_1\mathcal{K}\left(V\right)\text{d}\sigma }\right)-{\displaystyle {\int }_{\partial \Omega }{\left(g_1\right)}_t\mathcal{K}\left(V\right)\text{d}\sigma }$ (3.40)

Using the Trace Theorem (see for instance [23]) on the second term of the right hand side of (3.40), we obtain:

${\displaystyle {\int }_{\partial \Omega }{g}_1{\left(\mathcal{K}\left(V\right)\right)}_t\text{d}\sigma }\le \frac{\text{d}}{\text{d}t}\left({\displaystyle {\int }_{\partial \Omega }{g}_1\mathcal{K}\left(V\right)\text{d}\sigma }\right)+C{\Vert {\left(g_1\right)}_t\Vert }_{\infty }{\Vert \mathcal{K}\left(V\right)\Vert }_{H^1}$ (3.41)

Combining (3.36), (3.39) and (3.41), and using (2.24), we obtain

$\begin{array}{l}{\Vert \sqrt{D\left(V\right)}{V}_t\Vert }_{L^2}^2+\frac{1}{2}\frac{\text{d}}{\text{d}t}{\Vert \nabla \mathcal{K}\left(V\right)\Vert }_{L^2}^2\\ \le \frac{1}{2}{\Vert \sqrt{D\left(V\right)}{V}_t\Vert }_{L^2}^2+C{\Vert F\left(V\right)w_t\Vert }_{L^2}^2+C\left({\Vert w\Vert }_{L^{\infty }\left(L^{\infty }\right)}\right){\Vert \nabla \left(\mathcal{K}\left(V\right)\right)\Vert }_{L^2}^2\\ +\frac{\text{d}}{\text{d}t}\left(F\left(V\right)w,\nabla \left(\mathcal{K}\left(V\right)\right)\right)+\frac{\text{d}}{\text{d}t}\left({\displaystyle {\int }_{\partial \text{Ω}}{g}_1\mathcal{K}\left(V\right)\text{d}\sigma }\right)+C{\Vert {\left(g_1\right)}_t\Vert }_{\infty }{\Vert \mathcal{K}\left(V\right)\Vert }_{H^1}.\end{array}$ (3.42)

Next hide the first term on the right hand side of (3.42) in the like term on its left hand side, apply the Gronwall Lemma, and take the sup over the interval $\left[0,T_0\right]$ to get

$\begin{array}{l}{\Vert \sqrt{D\left(V\right)}{V}_t\Vert }_{L^2\left(L^2\right)}^2+{\Vert \nabla \left(\mathcal{K}\left(V\right)\right)\Vert }_{L^{\infty }\left(L^2\right)}^2\\ \le C\underset{0\le t\le T_0}{\sup }\left(\left|\left(F\left(V\right)w,\nabla \left(\mathcal{K}\left(V\right)\right)\right)\right|+{\Vert F\left(V\right)w_t\Vert }_{L^2\left(L^2\right)}^2\right)\\ +C\left(\left|\left({F\left(V^0\right)w|}_{t=0},\nabla \left(\mathcal{K}\left(V^0\right)\right)\right)\right|+{\Vert \nabla \left(\mathcal{K}\left(V^0\right)\right)\Vert }_{\left(L^2\right)}^2\right)\\ +\underset{0\le t\le T_0}{\sup }\left|{\displaystyle {\int }_{\partial \Omega }{g}_1\left(\mathcal{K}\left(V\right)\right)\text{d}\sigma }\right|+\left|{\displaystyle {\int }_{\partial \Omega }{g}_1\mathcal{K}\left(V^0\right)\text{d}\sigma }\right|\\ +C{\Vert {\left(g_1\right)}_t\Vert }_{\infty }{\displaystyle {\int }_0^{T_0}\left({\Vert \mathcal{K}\left(V\right)\Vert }_{H^1\left(\Omega \right)}\right)\left(\tau \right)\text{d}\tau },\end{array}$ (3.43)

Next, we bound the first term on the right hand side of (3.43).

$\begin{array}{l}C\underset{0\le t\le T_0}{\max }\left|\left(F\left(V\right)w,{\nabla \left(\mathcal{K}\left(V\right)\right)|}_{t=T_0}\right)\right|\\ \le \frac{1}{4}{\Vert \nabla \left(\mathcal{K}\left(V\right)\right)\Vert }_{L^{\infty }\left(L^2\right)}^2+C\left({\Vert w\Vert }_{\infty }\right){\Vert F\left(V\right)\Vert }_{L^{\infty }\left(L^2\right)}^2.\end{array}$ (3.44)

Using the Trace Theorem (see [23]), we handle the fifth term on the right hand side of (3.43) as follows.

(3.45)

where we have used the arithmetic-geometric inequality and chosen $ϵ>0$ so that

$C_1\frac{ϵ}{2}\le \frac{1}{4}$ .

With the help of Holder Inequality, the last term, on the right hand side of (3.43), can be dealt with as follows.

${\Vert {\left(g_1\right)}_t\Vert }_{\infty }{\displaystyle {\int }_0^{T_0}\left({\Vert \mathcal{K}\left(V\right)\Vert }_{H^1\left(\text{Ω}\right)}\right)\left(\tau \right)\text{d}\tau }\le C{\Vert \mathcal{K}\left(V\right)\Vert }_{L^2\left(H^1\right)}{\Vert g_t\Vert }_{\infty }{\Vert 1\Vert }_{L^2\left(H^1\right)}\le C\left({\Vert {\left(g_1\right)}_t\Vert }_{\infty }\right)$ (3.46)

by Lemma 3.5, if we assume that ${\left(g_1\right)}_t$ is bounded in $\partial \Omega \times \left(0,T_0\right]$ .

Also, by Lemma 3.5 and the regularity assumption on $F$ ,

${\Vert F\left(V\right)\Vert }_{L^{\infty }\left(L^2\right)}\le C,$ (3.47)

so, combining (3.43) through (3.47), and then hiding the appropriate terms, we get the lemma.

Remark 3.8

1. Unlike Lemma 3.5, Lemma 3.7 does use conditions (2.16) in its proof.

2. By (3.34), we have

$\begin{array}{c}{\Vert {\left(\mathcal{K}\left(V\right)\right)}_t\Vert }_{L^2\left(L^2\right)}={\Vert D\left(V\right)V_t\Vert }_{L^2\left(L^2\right)}={\Vert \sqrt{D\left(V\right)}\sqrt{D\left(V\right)}{V}_t\Vert }_{L^2\left(L^2\right)}\\ \le C{\Vert \sqrt{D\left(V\right)}{V}_t\Vert }_{L^2\left(L^2\right)}\le C.\end{array}$ (3.48)

Hence (3.33).

3. According to (3.30), $\mathcal{K}\left(V\right)$ is smoother than $V$ , the solution to (2.21), so this explains why, in numerical approximations of Problem 2.21, one often approximates first $\mathcal{K}\left(V\right)$ , and then uses the invertibility of the function $\mathcal{K}$ to find $V$ (see, for example [3] [4] [20] [24]).

We state and prove below another well-posedness theorem that gives uniqueness for $\mathcal{K}\left(V\right)$ in $L^{\infty }\left(0,T_0;L^2\left(\Omega \right)\right)$ and $L^2\left(0,T_0;H^1\left(\Omega \right)\right)$ , thus uniqueness for $V$ . A similar result was established for (1.4) in [2] for $\varphi \equiv 1$ .

Theorem 3.9 Under the conditions of Theorem 3.3

$\begin{array}{l}\underset{0\le t\le T_0}{\text{sup}}\left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right),V_2-V_1\right)+\eta {\Vert \nabla \left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\right)\Vert }_{L^2\left(L^2\right)}^2\\ \le C\left\{\left(\mathcal{K}\left(V_2^0\right)-\mathcal{K}\left(V_1^0\right),V_2^0-V_1^0\right)+{\Vert V_2^0-V_1^0\Vert }_{L^2}^{\frac{2}{2+\mu }}+{\Vert g_{1,2}-g_{1,1}\Vert }_{L^2\left(L^2\left(\partial \Omega \right)\right)}^2\right\}\end{array}$ (3.49)

with $\eta$ a positive constant, and $C=C\left({\Vert w\Vert }_{\infty },{\Vert w_t\Vert }_{L^{\infty }\left(L^2\right)}\right)$ .

Proof. In (3.11), replace $T\left(V_2-V_1\right)$ with $\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)$ , to get

$\begin{array}{l}\left(\frac{\partial \left(V_2-V_1\right)}{\partial t},\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\right)-\left(\left(F\left(V_2\right)-F\left(V_1\right)\right)w,\nabla \left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\right)\right)\\ +\left(\nabla \left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\right),\nabla \left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\right)\right)\\ ={\displaystyle {\int }_{\partial \Omega }\left(g_{1,2}-g_{1,1}\right)\left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\right)\text{d}\sigma },\end{array}$ (3.50)

which can be rewritten as

$\begin{array}{l}\left(\frac{\partial \left(V_2-V_1\right)}{\partial t},\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\right)+{\Vert \nabla \left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\right)\Vert }_{L^2}^2\\ =\left(\left(F\left(V_2\right)-F\left(V_1\right)\right)w,\nabla \left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\right)\right)\\ +{\displaystyle {\int }_{\partial \Omega }\left(g_{1,2}-g_{1,1}\right)\left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\right)\text{d}\sigma }.\end{array}$ (3.51)

Using the product rule, we can rewrite the first term on the left hand side of (3.51) as follows:

$\begin{array}{l}\left(\frac{\partial \left(V_2-V_1\right)}{\partial t},\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\right)\\ =\frac{\text{d}}{\text{d}t}\left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right),V_2-V_2\right)-\left({\left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\right)}_t,V_2-V_1\right).\end{array}$ (3.52)

Next, we rewrite (3.51) using (3.52).

$\begin{array}{l}\frac{\text{d}}{\text{d}t}\left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right),V_2-V_1\right)+{\Vert \nabla \left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\right)\Vert }_{L^2}^2\\ =\left(\left(F\left(V_2\right)-F\left(V_1\right)\right)w,\nabla \left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\right)\right)+\left({\left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\right)}_t,V_2-V_1\right)\\ +{\displaystyle {\int }_{\partial \Omega }\left(g_{1,2}-g_{1,1}\right)\left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\right)\text{d}\sigma },\end{array}$ (3.53)

or

$\begin{array}{l}\frac{\text{d}}{\text{d}t}\left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right),V_2-V_1\right)+{\Vert \nabla \left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\right)\Vert }_{L^2}^2\\ \le C{\Vert \left(F\left(V_2\right)-F\left(V_1\right)\right)w\Vert }_{L^2}^2+\frac{1}{4}{\Vert \nabla \left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\right)\Vert }_{L^2}^2\\ +{\Vert {\left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\right)}_t\Vert }_{L^{\gamma }}{\Vert V_2-V_1\Vert }_{L^{2+\mu }}+\frac{1}{2\epsilon }{\Vert g_{1,2}-g_{1,1}\Vert }_{L^2\left(\partial \text{Ω}\right)}^2\\ +C_2\frac{\epsilon }{2}{\Vert \mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\Vert }_{H^1}^2,\end{array}$ (3.54)

where $\gamma ,\epsilon >0$ , and

$\frac{1}{\gamma }+\frac{1}{2+\mu }=1$ (3.55)

with $\mu$ as in (2.16). Here we have used the Trace Theorem followed by the arithmetic-geometric inequality. Now

${\Vert \mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\Vert }_{H^1}^2={\Vert \mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\Vert }_{L^2}^2+{\Vert \nabla \left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\right)\Vert }_{L^2}^2,$ (3.56)

so choosing $\epsilon >0$ in (3.54) such that $C_2\frac{\epsilon }{2}=\frac{1}{4}$ , and then hiding the appropriate

terms from the left hand side to the righthand side, we obtain

$\begin{array}{l}\frac{\text{d}}{\text{d}t}\left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right),V_2-V_1\right)+\frac{1}{2}{\Vert \nabla \left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\right)\Vert }_{L^2}^2\\ \le C{\Vert \left(F\left(V_2\right)-F\left(V_1\right)\right)w\Vert }_{L^2}^2+{\Vert {\left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\right)}_t\Vert }_{L^{\gamma }}{\Vert V_2-V_1\Vert }_{L^{2+\mu }}\\ +C\left\{{\Vert g_{1,2}-g_{1,1}\Vert }_{L^2\left(\partial \Omega \right)}^2+{\Vert \mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\Vert }_{L^2}^2\right\}\end{array}$ (3.57)

after also hiding the second term on the righthand side of (3.54) in the like term on its left hand side.

Next, we integrate (3.57) over the interval $\left[0,t\right]$ , with $0<t\le T_0$ , and then take the supremum over $\left(0,T_0\right]$ to obtain:

$\begin{array}{l}\underset{0\le t\le T_0}{\text{sup}}\left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right),V_2-V_1\right)+\eta {\displaystyle {\int }_0^{T_0}{\Vert \nabla \left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\right)\Vert }_{L^2}^2\left(\tau \right)\text{d}\tau }\\ \le \frac{1}{2}{\displaystyle {\int }_0^{T_0}{\Vert \left(F\left(V_2\right)-F\left(V_1\right)\right)w\Vert }_{L^2}^2\left(\tau \right)\text{d}\tau }\\ +{\displaystyle {\int }_0^{T_0}{\Vert {\left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\right)}_t\Vert }_{L^{\gamma }}{\Vert V_2-V_1\Vert }_{L^{2+\mu }}\left(\tau \right)\text{d}\tau }\\ +C\left(\mathcal{K}\left(V_2^0\right)-\mathcal{K}\left(V_1^0\right),V_2^0-V_1^0\right)\\ +C\left\{{\Vert g_{1,2}-g_{1,1}\Vert }_{L^2\left(L^2\left(\partial \text{Ω}\right)\right)}^2+{\Vert \mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\Vert }_{L^2\left(L^2\right)}^2\right\}.\end{array}$ (3.58)

with $\eta$ a positive constant. By Holder inequality and (3.55),

$\begin{array}{l}{\displaystyle {\int }_0^{T_0}{\Vert {\left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\right)}_t\Vert }_{L^{\gamma }}{\Vert V_2-V_1\Vert }_{L^{2+\mu }}\left(\tau \right)\text{d}\tau }\\ \le {\Vert {\left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\right)}_t\Vert }_{L^{\gamma }\left(L^{\gamma }\right)}{\Vert V_2-V_1\Vert }_{L^{2+\mu }\left(L^{2+\mu }\right)}\end{array}$ (3.59)

Also, by (3.55), $\gamma \le 2$ , so that

${\Vert {\left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\right)}_t\Vert }_{L^{\gamma }\left(L^{\gamma }\right)}\le {\Vert {\left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\right)}_t\Vert }_{L^2\left(L^2\right)}\le C,$ (3.60)

by (3.33).

Finally, using (2.22), Theorem 3.3, and Corollary 3.4, we get Theorem 3.9.

Remark 3.10 Even if the porosity of the medium is not a function of $t$ , but still remains a function of the spacial variable, this analysis has the advantage of showing some explicit role of the porosity $\varphi$ in the advection (transport) through the

presence of the vector $w=\frac{\nabla \varphi }{\varphi }$ in the new equation.

4. The Full Saturation Equation

Up until now in this work, we have investigated the saturation equation (1.1) assuming that there is no transport term, i.e. the pure parabolic problem 1.7. With the change of the unknown function, a transport term was introduced and the transformed problem looked like a saturation problem (advection-diffusion problem). Now we consider the full saturation problem:

$\frac{\partial \left(\varphi S\right)}{\partial t}+\nabla \cdot \left(f\left(S\right)u\right)-\nabla \cdot \left(k\left(S\right)\nabla S\right)=Q\left(S\right)\text{in}\Omega \times \left(0,T_0\right)$ (4.1)

$\left(f\left(s\right)u-k\left(S\right)\nabla S\right)\cdot n=q_1\text{on}\partial \Omega \times \left(0,T_0\right)$ (4.2)

$S\left(x,0\right)=S_0\left(x\right)\text{in}\Omega ,$ (4.3)

where Ω is a sufficiently regular, bounded domain of $R^n$ , $n=1,2$ , $T_0>0$ , and $0\le S_0\left(x\right)\le 1$ on Ω. Here $u$ is the so-called Darcy velocity and $f$ is the

fractional flow function. We still assume that $\frac{\partial \varphi }{\partial t}$ is not necessarily 0.

One assumption often made on $f$ is

$f^{\prime }\left(0\right)=f^{\prime }\left(1\right)=0,$ (4.4)

see for instance [2]-[4].

4.1. The Transformed Problem

With the transformation, $V=\varphi S$ , made in section 2, if we put

$g\left(V\left(\cdot ,\cdot \right)\right)=f\left(\frac{V\left(\cdot ,\cdot \right)}{\varphi \left(\cdot ,\cdot \right)}\right)=f\left(S\left(\cdot ,\cdot \right)\right),$ (4.5)

equations (4.1) through (4.3) become

$\frac{\partial V}{\partial t}+\nabla \cdot \left(g\left(V\right)u+F\left(V\right)w\right)$

$-\Delta \mathcal{K}\left(V\right)=\tilde{Q}\left(V\right)\text{in}\Omega \times \left(0,T_0\right)$ (4.6)

$\left(g\left(V\right)u+F\left(V\right)w-D\left(V\right)\nabla V\right)\cdot n={\tilde{q}}_1\text{on}\partial \Omega \times \left(0,T_0\right)$ (4.7)

$V\left(x,0\right)=V_0\left(x\right),\text{in}\Omega ,$ (4.8)

where $V_0\left(x\right):=\varphi \left(x,0\right)S_0\left(x\right)$ , and $F$ , $D$ , and $\mathcal{K}$ are defined by (2.9), (2.8) and (2.19), respectively.

Recall that

$w=\frac{\nabla \varphi }{\varphi }.$ (4.9)

We make the more or less restrictive condition

$k^{\prime }\left(1\right)=0.$ (4.10)

However, this condition is not that restrictive, since, in the modeling of two-phase flow through porous media, conductivity functions that satisfy this condition are proposed (see for example [18]).

We define the vector $G$ by

$G\left(V\right)=g\left(V\right)u+F\left(V\right)w.$ (4.11)

Then, letting $Q\equiv 0$ , problem (4.6)~(4.8) becomes

$\frac{\partial V}{\partial t}+\nabla \cdot G\left(V\right)-\Delta \mathcal{K}\left(V\right)=0\text{in}\Omega \times \left(0,T_0\right)$ (4.12)

$\left(G\left(V\right)-D\left(V\right)\nabla V\right)\cdot n={\tilde{q}}_1\text{on}\partial \Omega \times \left(0,T_0\right)$ (4.13)

$V\left(x,0\right)=V_0\left(x\right),\text{in}\Omega .$ (4.14)

We have

$g^{\prime }\left(V\right)=\frac{1}{\varphi }{f^{\prime }\left(v\right)|}_{v=\frac{V}{\varphi }}.$ (4.15)

Therefore

$g^{\prime }\left(0\right)=\frac{1}{\varphi }{f}^{\prime }\left(0\right)=0,$ (4.16)

and

$g^{\prime }\left(\varphi \right)=\frac{1}{\varphi }{f}^{\prime }\left(1\right)=0,$ (4.17)

by (4.4). Thus (2.22) holds for $F$ replaced with $g$ and we have

${\left|g\left(v_2\right)-g\left(v_1\right)\right|}^2\le C^{\text{*}}\left(\mathcal{K}\left(v_2\right)-\mathcal{K}\left(v_1\right)\right)\left(v_2-v_1\right),$ (4.18)

where, for convenience, we call the constant $C^{\text{*}}$ again, quit to take the maximum of the two constants.

We note that if $w=0$ (that is the case when $\varphi \equiv 1$ ), then we get the saturation equation (1.4), and if $u=0$ , we get the pure parabolic equation (1.7).

We see that, for problem (4.12) - (4.14), conditions for Theorem 3.3 and Theorem 3.9 are satisfied, therefore the conclusions of these two theorems hold. We state these facts as corollaries.

Corollary 4.1 Under the conditions for Theorem 3.3, assume that (4.4) holds. Let $V_1$ and $V_2$ be solutions of Problem (4.12) - (4.14) corresponding to the initial conditions $V_1^0$ and $V_2^0$ respectively and the boundary conditions ${\tilde{g}}_{1,2}$ and ${\tilde{g}}_{1,2}$ respectively. Then

$\begin{array}{l}\underset{0\le t\le T_0}{\text{sup}}{\Vert V_2-V_1\Vert }_{{\left(H^1\right)}^{\text{*}}}^2+{\displaystyle {\int }_0^{T_0}\left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right),V_2-V_1\right)\left(\tau \right)\text{d}\tau }\\ \le C\left\{{\Vert V_2^0-V_1^0\Vert }_{{\left(H^1\right)}^{*}}^2+{\Vert {\tilde{g}}_{1,2}-{\tilde{g}}_{1,1}\Vert }_{L^2\left(L^2\left(\partial \Omega \right)\right)}^2\right\},\end{array}$ (4.19)

where $C=C\left(\varphi ,{\Vert u\Vert }_{\infty },{\Vert w\Vert }_{\infty }\right)$ .

Corollary 4.2 Under the conditions of Corollary 4.1, we have

$\begin{array}{l}\underset{0\le t\le T_0}{\sup }\left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right),V_2-V_1\right)+\eta {\Vert \nabla \left(\mathcal{K}\left(V_2\right)-\mathcal{K}\left(V_1\right)\right)\Vert }_{L^2\left(L^2\right)}^2\\ \le C\left\{\left(\mathcal{K}\left(V_2^0\right)-\mathcal{K}\left(V_1^0\right),V_2^0-V_1^0\right)+{\Vert V_2^0-V_1^0\Vert }_{L^2}^{\frac{2}{2+\mu }}+{\Vert {\tilde{g}}_{1,2}-{\tilde{g}}_{1,1}\Vert }_{L^2\left(L^2\left(\partial \Omega \right)\right)}^2\right\}\end{array}$ (4.20)

with $\eta$ a positive constant and $C=\left(\varphi ,{\Vert u\Vert }_{\infty },{\Vert u_t\Vert }_{L^2\left(L^{\infty }\right)},{\Vert w\Vert }_{\infty },{\Vert w_t\Vert }_{L^2\left(L^{\infty }\right)}\right)$ .

Also the conclusions of Lemma 3.5 and Lemma 3.7 hold for Problem (4.12) -(4.14) as well.

4.2. Back to the Initial Problem

Finally we go back to the initial saturation problem (4.1) - (4.3).

Because of (1.11) and (2.1), some of the estimates we proved for $V$ transfer directly to $S$ , the solution of (4.1) - (4.3). For example, we have:

${\Vert S_2-S_1\Vert }_{L^{\mu +2}\left(L^{\mu +2}\right)}^{\mu +2}={\Vert \frac{V_2}{\varphi }-\frac{V_1}{\varphi }\Vert }_{L^{\mu +2}\left(L^{\mu +2}\right)}^{\mu +2}\le {\left(\frac{1}{{\varphi }_0}\right)}^{\mu +2}{\Vert V_2-V_1\Vert }_{L^{\mu +2}\left(L^{\mu +2}\right)}^{\mu +2},$ (4.21)

by (1.11).

So that, by (3.22), we obtain the following.

Corollary 4.3 Suppose that $S_1$ and $S_2$ are solutions of Problem (4.1) - (4.3) corresponding to the initial conditions $S_1^0$ and $S_2^0$ , respectively. Assume $q_1\equiv 0$ . Then

${\Vert S_2-S_1\Vert }_{L^{\mu +2}\left(L^{\mu +2}\right)}^{\mu +2}\le C\left(\varphi ,{\Vert w\Vert }_{L^{\infty }\left(L^{\infty }\right)}\right){\Vert S_2^0-S_1^0\Vert }_{{\left(H^1\left(\Omega \right)\right)}^{\text{*}}}^2$ (4.22)

This resembles a result (established for (1.4) for $\varphi \equiv 1$ ) in [2], except that the constant depends now on $\varphi$ explicitly. making (4.22) more general than the one in [2].

Other estimates for $S$ will depend on the smoothness of $\varphi$ , through the vector $w$ and $w_t$ , for instance, when $\varphi$ is a function of $t$ or $x$ , or both.

5. Conclusion and Perspectives

In this paper, we have investigated the saturation problem (4.1) - (4.3), first by analyzing the case $u=0$ , through a change of unknown function, change that

introduced the term $w=\frac{\nabla \varphi }{\varphi }$ , which we mathematically interpret as a transport

term. A physical interpretation might come with a numerical treatment of the transformed problem. A priori estimates have been established as well as some regularity results. The latter results will be useful for regularization and numerical treatments of the problem. Because of the degenerate nature of the problem, we often opt for a regularization before a numerical treatment can be considered. In a forthcoming paper, we investigate a regularization of the problem in the line of [2] that follows the outline of a contribution chapter of the author in [19]. For the numerical approximation, a question naturally comes into mind: How will the introduction of a new transport term influence the performance of the various schemes proposed for the solution of the saturation problem (4.1) - (4.3)? Our perspective is to investigate this aspect of the problem in the very near future, as well as the full pressure/saturation system when the porosity is not independent of the time variable, based on error analyses for the proposed numerical methods that will follow in a sequel of this work.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Fadimba, K.B. (1993) Regularization and Numerical Methods for a Class of Porous Medium Equations. Ph.D. Thesis, University of South Carolina.
[2]	Fadimba, K.B. and Sharpley, R.C. (1995) A Priori Estimates and Regularization for a Class of Porous Medium Equations. The Nonlinear World, 2, 13-41.
[3]	Rose, M.E. (1980) Numerical Methods for Flows through Porous Media—II. Computers & Mathematics with Applications, 6, 99-122. https://doi.org/10.1016/0898-1221(80)90063-2
[4]	Rose, M.E. (1983) Numerical Methods for Flows through Porous Media. I. Mathematics of Computation, 40, 435-467. https://doi.org/10.1090/s0025-5718-1983-0689465-6
[5]	Smylie, D.L. (1989) A Near Optimal Order Approximation to a Class of Two-Sided Nonlinear Parabolic Partial Differential Equations. Ph.D. Thesis, University of Wyoming.
[6]	Vazquez, J.L. (2007) The Porous Medium Equation: Mathematical Theory. Oxford University Press.
[7]	Chavent, G. and Jaffre, J. (1986) Mathematical Models and Finite Element for Reservoir Simulation: Single Phase, Multiphase and Multicomponent Flows through Porous Media. North-Holland.
[8]	Bear, J. and Verruijt, A. (1987) Modeling Groundwater Flow and Pollution. D. Reidel Publication Company.
[9]	Ewing, R.E. (1983) 1. Problems Arising in the Modeling of Processes for Hydrocarbon Recovery. In: Ewing, R.E., Ed., The Mathematics of Reservoir Simulation, Society for Industrial and Applied Mathematics, 3-34. https://doi.org/10.1137/1.9781611971071.ch1
[10]	Fadimba, K.B. (2007) On Existence and Uniqueness for a Coupled System Modeling Immiscible Flow through a Porous Medium. Journal of Mathematical Analysis and Applications, 328, 1034-1056. https://doi.org/10.1016/j.jmaa.2006.06.012
[11]	Kieu, T. (2017) A Mixed Finite Element Approximation for Darcy-Forchheimer Flows of Slightly Compressible Fluids. Applied Numerical Mathematics, 120, 141-164. https://doi.org/10.1016/j.apnum.2017.05.006
[12]	Chen, Z., Huan, G. and Ma, Y. (2006). Computational Methods for Multiphase Flows in Porous Media. Society for Industrial and Applied Mathematics. https://doi.org/10.1137/1.9780898718942
[13]	Pop, I.S. (1998) Regularization Methods in the Numerical Analysis of Some Degnerate Parabolic Equations. Ph.D. Thesis, Babes-Bolyai University.
[14]	Arbogast, T. (1992) The Existence of Weak Solutions to Single Porosity and Simple Dual-Porosity Models of Two-Phase Incompressible Flow. Nonlinear Analysis: Theory, Methods & Applications, 19, 1009-1031. https://doi.org/10.1016/0362-546x(92)90121-t
[15]	Xiao, W. and Zhang, Y. (2022) Gevrey Regularity and Time Decay of Fractional Porous Medium Equation in Critical Besov Spaces. Journal of Applied Mathematics and Physics, 10, 91-111. https://doi.org/10.4236/jamp.2022.101008
[16]	Wojnar, R. (2010) Boussinesq Equation for Flow in an Aquifer with Time Dependent Porosity. Bulletin of the Polish Academy of Sciences: Technical Sciences, 58, 165-170. https://doi.org/10.2478/v10175-010-0017-3
[17]	Romero-Guzmán, E.T., Reyes-Gutiérrez, L.R., Klapp, J. and Hernández-Mendoza, H. (2018) Characterization of Natural Vuggy Fractured Porous Medium: Its Physicochemical Properties, Porosity and Permeability. Journal of Minerals and Materials Characterization and Engineering, 6, 38-49. https://doi.org/10.4236/jmmce.2018.61004
[18]	Fadimba, K.B. (2011) Pressure/Saturation System for Immiscible Two-Phase Flow: Uniqueness Revisited. Applied Mathematics, 2, 541-550. https://doi.org/10.4236/am.2011.25071
[19]	Fadimba, K.B. (2013) An Analysis for a Formulation of the Porous Medium Equation. In: Zhao, C.Y., Eds., Focus on Porous Medium Research, Material Science and Technologies, NOVA Science Publishers, 85-120.
[20]	Fadimba, K.B. and Sharpley, R.C. (2004) Galerkin Finite Element Method for a Class of Porous Medium Equations. Nonlinear Analysis: Real World Applications, 5, 355-387. https://doi.org/10.1016/j.nonrwa.2003.07.001
[21]	Ciarlet, P.G. and Oden, J.T. (1978) The Finite Element Method for Elliptic Problems. Journal of Applied Mechanics, 45, 968-969. https://doi.org/10.1115/1.3424474
[22]	Glowinski, R. (1984) Numerical Methods for Nonlinear Variational Problems. Springer-Verlag.
[23]	Brenner, S.C. and Scott, L.R. (1994) The Mathematical Theory of Finite Element Methods. Springer.
[24]	Fadimba, K.B. (2005) A Linearization of a Backward Euler Scheme for a Class of Degenerate Nonlinear Advection-Diffusion Equations. Nonlinear Analysis: Theory, Methods & Applications, 63, e1097-e1106. https://doi.org/10.1016/j.na.2004.09.038