Electronic Journal of Differential Equations, Vol. 2022 (2022), No. 06, pp. 1–22.

ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu

MATHEMATICAL ANALYSIS OF A

DUPUIT-RICHARDS MODEL

SAFAA AL NAZER, CAROLE ROSIER, MUNKHGEREL TSEGMID

Abstract. This article concerns an alternative model to the 3D-Richards

equation to describe the flow of water in shallow aquifers. The model couples
the two dominant types of flow existing in the aquifer. The first is described

by the classic Richards problem in the upper capillary fringe. The second

results from Dupuit’s approximation after vertical integration of the conser-
vation laws between the bottom of the aquifer and the saturation interface.

The final model consists of a strongly coupled system of parabolic-type partial

differential equations that are defined in a time-dependent domain. First, we
show how taking the low compressibility of the fluid into account eliminates

the nonlinearity in the time derivative of the Richards equation. Then, the

general framework of parabolic equations is used in non-cylindrical domains
to give a global in time existence result to this problem.

1. Introduction

Populated areas are increasingly affected by contamination of soil and ground-
water. Many modeling have been developed to study the vulnerability of aquifers
to agricultural pollution, with a particular focus on the supply of nitrates. There
is a wide variety of processes involved (chemical, hydrogeological, anthropic, . . . )
acting in a wide range of temporal and geometrical length scales. But we notice
that the main point for the derivation of the hydrogeological model is linked to a
good description of the flow between the ground level (the level of the anthropic
processes) and the water table. This will be crucial when studying the transport
of chemical components in the aquifer. It turns out that many chemical reactions
are expected in the first meters of the subsoil, where oxygen is still very present.
In particular, chemical species that reach the water table are not necessarily the
same as those that have left the surface. This yields different speeds of the reactive
kinetics. As a result, for an efficient mathematical modeling, the time upscaling
process in this zone must keep track of all the time scales.

Only the hydrogeological question will be considered in this work. Aquifers are
often characterized by a form of stratification of flows which enables the definition
of interfaces. The slowness of the natural dynamics ensures a smooth and stable
behavior for the interfaces. Besides, due to the dimensions of the aquifer, the flow
can be assumed essentially orthogonal to the equipotential (Dupuit’s hypothesis).

2020 Mathematics Subject Classification. 35A01, 35D30, 35K59, 35Q86, 35R35.
Key words and phrases. Dupuit-Richards equations; free boundary problems; global solution;

weak solution; fluid flow modeling.
©2022. This work is licensed under a CC BY 4.0 license.

Submitted May 5, 2021. Published January 17, 2022.

1



2 S. AL NAZER, C. ROSIER, M. TSEGMID EJDE-2022/06

The vertical integration of the Richards equation is thus possible, at least in the
saturated zone. In this spirit, many 2D models have been developed and used since
the 1960s (see for example the works of Jacob Bear, [5, 6]). For more historical notes
on the origin of groundwater modeling, we refer interested readers to [14, 15, 17, 24].
But the approach by vertical integration is only valuable for very precise length and
time scales, the time scale in particular being completely different from the typical
durations of chemical reactions. However, such 2D models are widely used, although
it is particularly difficult to correctly couple them to the flow in the unsaturated part
of the basement. Several numerical studies have been conducted in this direction.
Let us mention the work of [20] where the integrated model is directly coupled with
a surface model. In [7] and [30], the coupling of the surface and underground flows is
done using a Richards equation associated with a Signorini boundary condition (for
the surface behavior). A class of models is proposed in [8], which consists in coupling
purely vertical models (to describe the flow at a small time scale) with an horizontal
model (describing the flow at a long time scale). They admit the same behavior
than the 3D-Richards model for any time scale when the aquifer presents a small
deepness compared to its large horizontal dimensions. They describe the essentially
horizontal flow of a water table and the essentially vertical water supply flux from
the surface through the unsaturated part between the groundwater and the ground
level. In [27] a presentation of a rather similar model can be found, coupling 1D-
Richards equation with a simplified model in the saturated part. Finally, in [1], this
kind of model is integrated into a computational code called ”SHE” (for ”European
Hydrological System” and later became SHETRAN), in the case where the water
table remains away from ground level.

In this article, a model belonging to the “Dupuit-Richards” model class is ana-
lyzed. Indeed, the 3D-Richards equation is considered in the capillary fringe while
a vertical average of the mass conservation law is made in the saturated zone of
the aquifer. Pressure and normal fluxes transmission properties are imposed at the
saturation interface.

This model differs slightly from the one described in [8], firstly because the
complete 3D Richards equations are considered in the unsaturated part and not
only the vertical component of the flow. Nevertheless, the main difference lays in
the consideration of the low compressibility of the fluid. This makes it possible to
perform a change of variable which eliminates the nonlinearity of the time derivative
of the moisture content, but also to treat the degeneracy in the equation governing
the horizontal flow in the saturated part of the aquifer. On the other hand, the
coupling of flows between the two areas of the aquifer results from the continuity
property of the normal component of the flux at the interface of saturation, ensuring
the mass conservation at the interface.

Of course, the numerical behavior of this model should be similar to the one
obtained in [8], especially when the horizontal component of the hydraulic con-
ductivity tends towards zero in the unsaturated part of the aquifer. Finally, the
coupling of the 1D Richards equation with the 2D Dupuit approximation numer-
ically justifies this model since it is less expensive than the complete resolution of
the 3D Richards equation in terms of CPU time.

The mathematical study of the model is particularly delicate because of non-
linearities, the free boundary between each area of the aquifer and the difficulty
resulting from the coupling between the two zones, which is expressed here in terms



EJDE-2022/06 DUPUIT-RICHARDS MODEL 3

of flux at the interface. Finally, there is a general mathematical complexity in the
structure of the set of PDEs modeling the dynamics of underground water. Indeed,
when considering a free water table, we must face the gradual disappearance of
water in the desaturation zone and thus the disappearance of a main unknown of
the problem hence the mathematical challenges inherent in Richards equations.

There exists a huge body of literature regarding the classical Richards equations.
Let us mention the works of Alt et al ([3, 4]) and the papers [10, 18, 28] devoted
to the study of the degenerate in time equation

∂tθ(p)−∆p = 0,

where θ(p) denotes the moisture content. In the one-dimensional case, we also
quote the work of Yin ([33]) concerning the existence of weak solution for the fully
degenerate problem

∂tθ(p)− ∂x(κ(θ(p))∂xp) = 0,

when just assuming that θ′, κ′ > 0.
Classically, the Kirchoff transform is applied to the Richards equation (under

appropriate assumptions about porosity and permeability) to eliminate the nonlin-
earity in the diffusive term. In this work, the hypothesis of low compressibility of
water is exploited instead, to eliminate the nonlinearity in the time derivative of the
Richards equation. Even if this term is assumed to be very low, its very existence
allows to define the one-to-one transformation that absorbs the degeneracy of the
moisture content. This transformation brings us back to the framework of quasilin-
ear parabolic equations on non-cylindrical domains to which the auxiliary domain
method introduced by Lions and Mignot [22, 25] can be applied to deal with the
free boundary. Regarding the flow in the saturation zone, the vertical mean of con-
servation laws in this part leads to a degenerate elliptic equation whose degeneracy
depends on the thickness of the zone. In addition, taking the compressibility of
water into account introduces degeneracy (in the time derivative) also depending
on the thickness of the saturated zone. A change of variable then allows to absorb
the two degenerated terms and return to a regular parabolic equation.

This article is organized as follows: In section 2, the model coupling 3D Richards
equation with the Dupuit horizontal approximation is introduced; consequences
taking the compressibility of the fluid into account in the modeling are particularly
detailed. Then a global in time existence result is given in Section 3 as well as
preliminary results about the auxiliary domains method. The proof of the Theorem
is performed in section 4. It consists in a fixed point strategy to deal with the
difficulties associated with nonlinearities and coupling.

2. Derivation of the model

The basis of the modeling is the mass conservation law written for fresh water
coupled with the classical Darcy law for porous media. First, the Richards equations
are obtained by taking the water compressibility into account. Then, the case of the
unsaturated zone is distinguished from that of the saturated zone as explained in
the introduction. For the three-dimensional description, we denote by x := (x, z),
x = (x1, x2) ∈ R2, z ∈ R, the usual coordinates.

2.1. Geometry. The aquifer is represented by a three-dimensional domain Ω :=
Ωx × (hbot, hsoil), Ωx ⊂ Rn with n ≥ 2 (x = (x1, x2)), function hbot (respect. hsoil)
describing its lower (respect. upper) topography. The upper and lower surfaces



4 S. AL NAZER, C. ROSIER, M. TSEGMID EJDE-2022/06

are thus defined by the graph of the functions hbot = hbot(x) and hsoil = hsoil(x),
x ∈ Ωx. We assume that

hsoil(x) > hbot(x), ∀x ∈ Ωx. (2.1)

More precisely the domain is given by

Ω =
{

(x, z) ∈ Ωx × R : z ∈]hbot(x), hsoil(x)[
}
. (2.2)

We always denote by ~ν the outward unit normal and ~e3 is the unitary vertical vector
pointing up. The boundary ∂Ω of Ω is divided into three zones (bottom, top and
vertical)

∂Ω = Γbot t Γsoil t Γver ,

with

Γbot :=
{

(x, z) ∈ Ω : z = hbot(x)
}
, Γsoil :=

{
(x, z) ∈ Ω : z = hsoil(x)

}
,

Γver :=
{

(x, z) ∈ Ω : x ∈ ∂Ωx
}
.

The model divides the flow description into two subregions of Ω (possibly time-
dependent) in each of which the flow presents a different behavior. We denote by h
the depth of the free interface separating the freshwater layer and the unsaturated
part of the aquifer. The definition of these zones is thus based on the function
h = h(t, x) which is an unknown of our problem. Then, for a given function
h = h(t, x) such that hbot ≤ h ≤ hsoil, we introduce:

Ω−t :=
{

(x, z) ∈ Ω : z < h(t, x)
}
, Ωt :=

{
(x, z) ∈ Ω : z > h(t, x)

}
, (2.3)

Γt :=
{

(x, z) ∈ Ω : z = h(t, x)
}
. (2.4)

2.2. Conservation laws. In view of the (large) dimensions of an aquifer (related
to the characteristic size of the porous structure of the underground), we consider
a continuous description of the porous medium. The effective velocity q of the flow
is thus related to the pressure P through the Darcy law associated with a nonlinear
anisotropic conductivity

q = −κ(P )K0

µ
(∇P + ρg∇z),

where ρ and µ are respectively the density and the viscosity of the fluid, K0 is the
permeability of the soil, κ(P ) is the relative conductivity and g the gravitational
acceleration constant. Introducing the hydraulic head H defined by

H =
P

ρ0g
+ z, (2.5)

the previous equation is written as follows,

q = −K∇H − κ(P )K0

µ
(ρ− ρ0)g∇z, K =

κ(P )K0ρ0g

µ
. (2.6)

In this relation, the matrix K is the hydraulic conductivity which expresses the
ability of the underground to conduct the fluid. ρ0 denotes the reference density
of the fluid. Next, the conservation of mass during displacement is given by the
equation

∂t(θρ) +∇ · (ρq) = ρQ, (2.7)



EJDE-2022/06 DUPUIT-RICHARDS MODEL 5

where Q denotes a generic source term (for production and replenishment). The
function θ is the volumetric moisture content defined by

θ = φs,

where φ is the porosity of the medium and s is the saturation. If the air present in
the unsaturated zone is assumed to have infinite mobility, the saturation s, and then
the function θ are thus considered monotone functions depending on the pressure
as we will detail it latter.

2.3. State equation for the fluid compressibility. The fluid is considered com-
pressible by assuming that pressure P is related to the density ρ as follows (cf. [11]):

dρ

ρ
= αP dP ⇔ ρ = ρ0e

αP (P−P0). (2.8)

The real number αP ≥ 0 represents the fluid compressibility coefficient and P0 is the
pressure of reference. Further assuming αP = 0 we would recover the incompressible
case. This compressibility coefficient αP is of course very low, but, as mentioned
in the introduction, the fact of not completely neglecting it makes it possible to
facilitate the mathematical analysis of the Richards equations.

2.4. Permeability tensor K0. The nonlinear hydraulic conductivity K is given

by K = κ(P ) ρ0g
µ K0. The soil transmission properties are characterized by the

porosity function φ and the permeability tensor K0(x, z). The matrix K0 is a 3× 3
symmetric positive definite tensor which describes the conductivity of the saturated
soil at the position (x, z) ∈ Ω. We introduce Kxx ∈ M22(R), Kzz ∈ R∗ and
Kxz ∈M21(R) such that

K0 =

(
Kxx Kxz

KT
xz Kzz

)
. (2.9)

2.5. Hypothesis. Let us now list the assumptions on the fluid and medium char-
acteristics but also on the flow which are meaningful in the context of this problem.

2.5.1. Hypothesis on the fluid and on the medium. In the model, the effects of the
rock compressibility are neglected, the porosity of the medium φ do not depend on
the pressure variations and it is thus assumed to be a constant.
Compressibility of the fluid. The fluid (namely here fresh water) is weakly
compressible. It means that

αP � 1. (2.10)

Let us exploit this assumption. In natural conditions and especially in an aquifer,
one observes small fluid mobility (defined by the ratio κ/µ). First consequence of
the low compressibility of the fluid combined with the low mobility of fluid appears
in the momentum equation. We perform a Taylor expansion with regard to P of the
density ρ in the gravity term of the Darcy equation. Neglecting the terms weighted
by αPκ/µ� 1 in (2.6), we obtain

q = −K∇H, K =
κ(P )ρ0g

µ
K0. (2.11)

Second consequence is ∇ρ · q � 1 which leads to the following simplification in the
mass conservation equation (2.7),

ρ∂tθ + θ∂tρ+ ρ∇ · q = ρQ.



6 S. AL NAZER, C. ROSIER, M. TSEGMID EJDE-2022/06

Neglecting in this way the variation of density in the direction of flow is sometimes
considered as an extra assumption called Bear’s hypothesis (cf [2]). Including (2.8),
that is ∂tρ = ραP∂tP in the latter equation, we obtain

ρ∂tθ + ρθαP∂tP + ρ∇ · q = ρQ.

After simplification by ρ > 0, we finally obtain

∂tθ + θαP∂tP +∇ · q = Q. (2.12)

Equivalently, using the hydraulic head (2.5) and the Darcy law (2.11), (2.12) can
be written

∂tθ + S0∂tH −∇ · (K∇H) = Q where S0 = ρ0gφαP . (2.13)

We notice that if the fluid is assumed incompressible, αP = 0, then (2.12) is the
classical Richards equation in pressure formulation. An adequate definition of the
volumetric moisture content θ and of the mobility function κ is the key of the model.
Richards hypothesis. The Richards model is moreover based on the assumption
that the air pressure in the underground equals the atmospheric pressure, thus is not
an unknown of the problem. One thus assumes that the saturation and the relative
conductivity of the soil are given as functions of the fluid pressure P , denoted
respectively by s = s(P ) and κ = κ(P ). We introduce the saturation pressure Ps
which is a fixed real number. The fully-saturated part of the medium corresponds
to the region {x;P (·,x) > Ps}, while it is partially-saturated in the capillary fringe
{x;Pd < P (·,x) ≤ Ps}. The dry part is defined by the set {x;P (·,x) ≤ Pd}. The
moisture content is such that

φ if P (·,x) > Ps (saturated zone),

θ(P ) if Pd < P (·,x) ≤ Ps with 0 ≤ θ(P ) ≤ φ and θ′(P ) > 0),

θ0 = φs0 if P (·,x) ≤ Pd (dry zone),

(2.14)

where s0 > 0 corresponds to a residual saturation which is positive. The associated
relative hydraulic mobility is then defined by

κ(P ) =


1 if P (·,x) > Ps (saturated zone),

κ(θ(P )) if Pd < P (·,x) ≤ Ps with 0 ≤ κ(P ) ≤ 1 and (κoθ)′(P ) > 0),

0 if P (·,x) ≤ Pd (dry zone).

(2.15)
There is a large choice of available models for s and κ. The most classical examples
for an air-water system are the van Genuchten model [32] with no-explicit depen-
dance on the bubbling pressure but with fitting parameters, and the Brooks and
Corey model [9]. The important point is that these models are such that

s(P ) = 1 ⇐⇒ P ≥ Ps,
κ(P ) = 1 ⇐⇒ P ≥ Ps.

(2.16)

In particular, the water pressure is greater than the bubbling pressure Ps if and
only if the soil is completely saturated.

2.5.2. Hypothesis on the flow. The following assumption is introduced for upscaling
the 3D problem to a 2D model in the saturated part of the domain.
Dupuit approximation (hydrostatic approach) Dupuit assumption consists
in considering that the hydraulic head is constant along each vertical direction
(vertical equipotentials). It is legitimate since one actually observes quasi-horizontal



EJDE-2022/06 DUPUIT-RICHARDS MODEL 7

displacements when the thickness of the aquifer is small compared to its width and
its length and when the flow is far from sinks and wells.

2.6. Model coupling vertical 3d-Richards flow and Dupuit horizontal
flow. From the compressible Richards equation (2.12), we will deduce the equations
governing the flow in each zone of the aquifer.
• Three-dimensional Richards equation in the upper capillary fringe. In
the unsaturated part of the aquifer, Ωt, the 3D Richards equation (2.12) holds

∂tθ + θαP∂tP +∇ · q = Q for (t, x, z) ∈ (0, T )× Ωt,

q · ~ν = 0 for (t, x, z) ∈ (0, T )×
(
Γrmsoil ∪ Γver

)
,

P
(
t, x, h(t, x)

)
= Ps for (t, x) ∈ (0, T )× Ωx,

P (0, x, z) = Pinit(x, z) for (x, z) ∈ Ω0.

(2.17)

The effective velocity q is given by

q = −K∇(
P

ρog
+ z), K =

κ(P )K0ρ0g

µ
.

We emphasize that the model (2.17) depends by definition on the depth h which is
expected to belong to the interval (hbot, hsoil).
• Dupuit horizontal flow in the saturated zone. In the saturated part of the
aquifer, Ω−t , the vertical average of the 3D Richards equation (2.13) describes the
horizontal flow of this part, thus reducing the 3D problem to a 2D problem. For
the upscaling procedure, we proceed as it was done in the context of the seawater
intrusion in [11, 12].

Upscaling procedure The vertical integration is performed between depths hbot and
h. Since θ(P ) = φ in the saturated zone, the vertical average (2.13) leads to∫ h

hbot

(
S0∂tH +∇ · q

)
dz =

∫ h

hbot

Qdz.

Bf = h − hbot denotes the thickness of the saturated zone and Q̃ the source term
representing distributed surface supply of fresh water into the free aquifer:

Q̃ =
1

Bf

∫ h

hbot

Qdz.

Applying Leibnitz rule to the first term in the left-hand side yields:∫ h

hbot

S0∂tHdz = S0
∂

∂t

∫ h

hbot

Hdz − S0H|z=h∂th+ S0H|z=hbot
∂thbot.

We denote by H̃ the vertically averaged hydraulic head

H̃ =
1

Bf

∫ h

hbot

Hdz.

Because of Dupuit approximation, H(x1, x2, z) ' H̃(x1, x2), x = (x1, x2) ∈ Ω,
z ∈ (hbot, h), hence we have∫ h

hbot

S0∂tHdz = S0Bf∂tH̃.



8 S. AL NAZER, C. ROSIER, M. TSEGMID EJDE-2022/06

Also ∫ h

hbot

∇ · q dz = ∇′ · (Bf q̃′) + q|z=h− · ∇(z − h)− q|z=h+
bot
· ∇(z − hbot),

where ∇′ = (∂x1
, ∂x2

), q′ = (qx1
, qx2

).

The averaged Darcy velocity q̃′ = 1
Bf

∫ h
hbot

q′ dz is

q̃′ = − 1

Bf

∫ h

hbot

(K∇′H) dz = −K̃∇′H̃, K̃ =
1

Bf

∫ h

hbot

K0ρ0g

µ
dz,

(we remind that κ(P ) = 1 for z ∈ (hbot, h)). The averaged mass conservation law
for the freshwater in the saturated zone finally reads

S0Bf∂tH̃ = ∇′ · (Bf K̃∇′H̃) + q|z=h+
bot
· ∇(z − hbot)

− q|z=h− · ∇(z − h) +Bf Q̃.
(2.18)

In this equation, the term Bf K̃ may be viewed as the dynamic transmissivity of
the freshwater layer. At this stage, we have obtained an undetermined system of
two PDEs (2.17)-(2.18) with three unknowns P , H̃ and h.
Fluxes and continuity equations across the interface

Our aim is now to include in the model the continuity and transfert properties
across interfaces. As a consequence, we express the two flux terms appearing in
(2.18) and the number of unknowns is reduced.
. Flux across the saturation interface: The saturation interface is characterized by
the cartesian equation

F (x1, x2, z, t) = 0 ⇐⇒ z − h(x1, x2, t) = 0,

so the unit normal vector ~ν to the interface is colinear to ∇(z − h). The relation
ruling continuity of the normal component of the velocity thus reads(

q|z=h+ − q|z=h−
)
· ~ν = 0 ⇐⇒ q|z=h+ · ∇(z − h) = q|z=h− · ∇(z − h). (2.19)

. Approximation of the flux q|z=h+ · ∇(z − h): The flux q|z=h+ · ∇(z − h) ex-
presses mass transfers between the two parts of the aquifer. As it is done in [8], we
approximate the flux by

q|z=h+ · ∇(z − h) '
∫ hsoil(x)

h(t,x)

(
φ
∂s(P )

∂t
+ φs(P )αP

∂p

∂t
−Q

)
dz. (2.20)

This approximation comes from the hypothesis of an almost null horizontal hy-
draulic conductivity (i.e. Kxx ' (0)) in the capillary fringe. This corresponds to
a flow almost vertical in this part of the aquifer. So the 3D-Richards equation is
reduced to an 1D-equation. Integrating this 1D equation between h and hsoil leads
to the approximation (2.20).

It is an essential difference with the mathematical analysis presented in [31]
in which the exchanges between the two parts of the aquifer were simplified and
represented by the addition of an external source term, thus decoupling the two
problems.
. Impermeable layer at z = hsoil: Since the lower layer is impermeable, there is no
flux across the boundary z = hbot:

q(hbot) · ∇(z − hbot) = 0. (2.21)



EJDE-2022/06 DUPUIT-RICHARDS MODEL 9

. Continuity equations: The continuity relation imposed on the interface enables
to properly reduce the number of unknowns in equations (2.17)-(2.18). Dupuit

approximation reads H̃ ' H|z=h− , the pressure P thus satisfies in Ω−t

P (t, x, z) = ρ0g
(
H̃(t, x)− z

)
for t ∈ [0, T [ (x, z) ∈ Ω−t . (2.22)

Also, the pressure is continuous across Γt, it follows that

P (t, x, h−) = P (t, x, h+) = Ps ⇐⇒ H̃ =
Ps
ρ0g

+ h. (2.23)

Equation (2.23) allows to substitute H̃ by h in Eq. (2.18), so we have

S0Bf∂th−∇′ · (Bf K̃∇′h)

= Bf Q̃−
∫ hsoil(x)

h(t,x)

(
φ
∂s(P )

∂t
+ φs(P )αP

∂p

∂t
−Q

)
dz in (0, T )× Ωx,

(2.24)

K̃∇′h · ~ν = 0 on (0, T )× ∂Ωx, (2.25)

with

Bf = (h− hbot), K̃ =
1

Bf

∫ h

hbot

K0ρ0g

µ
dz, S0 = ρ0gφαP . (2.26)

The homogeneous Neumann condition on ∂Ωx is assumed to simplify the presenta-
tion.

The final model (M) coupling 3D-Richards flow and Dupuit horizontal flow is
completed with initial and boundaries conditions. It thus consists in the following
system:
• In Ωt the 3d-Richards equation holds,

∂tθ(P ) + θαP∂tP +∇ · q = Q in (0, T )× Ωt,

q · ~ν = 0 on (0, T )×
(
Γrmsoil ∪ Γver

)
,

P
(
t, x, h(t, x)

)
= Ps in (0, T )× Ωx,

P (0, x, z) = P0(x, z) in Ω0,

where the saturation pressure Ps is assumed to be constant with respect to the time
and to the space. The effective velocity q is

q = −K∇(
P

ρog
+ z), K =

κ(P )K0ρ0g

µ
.

• In Ω−t the pressure P satisfies

P (t, x, z) = ρ0g
( Ps
ρ0g

+ h− z
)

in (0, T )× Ω−t .

• The depth of Γt, h in Ωx, satisfies

S0Bf∂th−∇′ · (Bf K̃∇′h) = Bf Q̃−
∫ hsoil(x)

h(t,x)

(
φ
s(P )

∂t
+ φ s(P )αP

∂p

∂t
−Q

)
dz,

K̃∇′h · ~ν = 0 on (0, T )× ∂Ωx,

h(0, x) = h0(x) in Ωx.



10 S. AL NAZER, C. ROSIER, M. TSEGMID EJDE-2022/06

3. Mathematical setting and main results

Problem (M) being a problem with free boundary, we are going to define the
general framework of parabolic equation in non cylindrical domains, introduced by
Lions and Mignot respectively in [22] and [25]. There are obviously many other
techniques to deal with free boundary problems, we refer interested readers to the
recent paper [19] and references therein.

3.1. Notation and auxiliary results. For any T > 0, let OT be the open domain
of R+ × Ω defined by

OT =
{

(t, x, z) ∈ (0, T )× Ω : h(t, x) < z
}

where h is the position of the interface Γt. We set

Ωt =
{

(x, z) ∈ Ω : z ∈]h(t, x), hsoil[
}
, OcT =

(
(0, T )× Ω

)
\ OT ,

Γ = ∂OT (boundary of OT ), Γ′ = Γ \ (Ω0 ∪ ΩT ) (lateral boundary of OT ).

We define

H0,1(OT ) = {u ∈ L2(OT ) : Dp
xu ∈ L2(OT ) for |p | ≤ 1},

where

Dp
xu = {Dα

xu : α = (α1, α2, α2) with |α| = p}.
It is an Hilbert space endowed with the norm

‖u‖H0,1(OT ) =
(

Σp≤1

∫
OT
|Dpu| dx dt

)1/2

.

F (OT ) denotes the closure in H0,1(OT ) of functions of D(ŌT ) null in a neighbor-
hood of Γt and F ′(OT ) its topological dual. Besides, we introduce

B(OT ) = {u ∈ F (OT ) :
du

dt
∈ F ′(OT )},

endowed with the Hilbertian norm

‖ · ‖2B(OT ) = ‖ · ‖2F (OT ) + ‖∂t · ‖2F ′(OT ).

Finally, B0(OT ) (resp. BT (OT )) is the closure in B(OT ) of functions of B(OT )
null in a neighborhood of t = 0 (resp. t=T). We now state some auxiliary results
proved in [22]

Lemma 3.1. If OT is sufficiently regular, we have

1. H0,1(OT ) = L2([0, T ];H1(Ωt)) where

L2(0, T ;H1(Ωt)) = {u(t, ·) ∈ H1(Ωt), t ∈ [0, T ], a.e. and ‖u‖H0,1(OT ) < +∞},

with ‖u‖H0,1(OT ) =
∫ T

0
‖u‖2H1(Ωt)

dt.

We have a similar result holds for F (OT ).

2. For u ∈ F (OT ), we can define γ(u), the trace of u on Γ′ in L2(Γ′).

Moreover u ∈ F (OT ) ⇐⇒ γ(u) = 0 on Γt.

3. Let u ∈ B(OT ), thus u ∈ BT (OT ) ⇐⇒ u(T, .) = 0.
4. for all u, v ∈ B(Os), we have

〈∂u
∂t
, v〉F ′,F + 〈∂v

∂t
, u〉F ′,F = (u(s, ·), v(s, ·))L2(Ωs) − (u(0, .), v(0, .))L2(Ω0). (3.1)



EJDE-2022/06 DUPUIT-RICHARDS MODEL 11

For the sake of brevity, we shall write H1(Ω) = W 1,2(Ω) and

V (Ω) = H1
0,Γbot

(Ω) = {u ∈ H1(Ω), u = 0 on Γbot}, V ′(Ω) = (H1
0,Γbot

(Ω))′.

The embeddings V (Ω) ⊂ L2(Ω) ⊂ V ′(Ω) are dense and compact. For any T > 0,
let W0(0, T,Ω) denotes the space

W0(0, T,Ω) :=
{
ω ∈ L2(0, T ;V (Ω)), ∂tω ∈ L2(0, T ;V ′(Ω))

}
,

endowed with the Hilbertian norm ‖·‖2W0(0,T,Ω) = ‖·‖2L2(0,T ;V (Ω))+‖∂t·‖
2
L2(0,T ;V ′(Ω)).

The following embeddings are continuous [23, prop. 2.1 and thm 3.1, chapter 1]

W0(0, T,Ω) ⊂ C([0, T ]; [V (Ω), V ′(Ω)] 1
2
) = C([0, T ];L2(Ω))

while the embedding

W (0, T,Ω) ⊂ L2(0, T ;L2(Ω)) (3.2)

is compact (Aubin’s Lemma, see [29]).
In the same way, we introduce the space

W (0, T,Ωx) :=
{
ω ∈ L2(0, T ;H1(Ωx)) : ∂tω ∈ L2(0, T ; (H1(Ωx))′)

}
,

endowed with the Hilbertian norm

‖ · ‖2W (0,T,Ωx) = ‖ · ‖2L2(0,T ;H1(Ωx)) + ‖∂t · ‖2L2(0,T ;(H1(Ωx)))′ .

The same compacity results hold true in this case.

3.2. Main results. We aim giving an existence result of physically admissible
weak solutions for model (M) completed by initial and boundary conditions. Let
us first detail the mathematical assumptions. We begin with the characteristics of
the porous structure. The study is limited to the isotropic case so K0 is assumed
to be a scalar. In the saturated part, the averaged hydraulic conductivity K̃ is
thus equal to the constant K0ρ0g/µ. Without loss of generality, we will assume a
zero source term, Q = 0. The initial data P0 ∈ H2(Ω) satisfies the compatibility
condition

P0(x, h0) = Ps in Ω0.

Let δ ∈ R be a positive number, we assume that h0 ∈ L∞(Ωx) is such that

hbot + δ ≤ h0 ≤ hsoil a.e. in Ωx. (3.3)

Functions θ and κ are pressure-dependent and we assume

θ ∈ C1(R), 0 < θ− := φs0 ≤ θ(x) ≤ θ+, θ′(x) ≥ 0 ∀x ∈ R, (3.4)

κ ∈ C(R), 0 < κ− ≤ κ(x) ≤ κ+ ∀x ∈ R. (3.5)

Before stating the main result of this work, we will transform the original problem
and bring us back to the framework introduced in [25].

The above assumptions on the fluid and the medium allow to eliminate the
nonlinearity in time of (2.17), namely assumptions (3.4)-(3.5) are sufficient to define
the primitive function P such that

P(P ) = θ(P ) + αP

∫ P

θ(s) ds.

A direct computation gives P ′(P ) = θ′(P ) + αP θ(P ) > αP θ− > 0, indeed by
previous hypothesis, we have θ′(P ) ≥ 0 and θ(P ) > φs0. Since θ ∈ C1(R), there



12 S. AL NAZER, C. ROSIER, M. TSEGMID EJDE-2022/06

exists θ′+ > 0 such that 0 ≤ θ′ ≤ θ′+ on the interval [Pd, Ps]. Since P is a bijective
application, the existence of p such that

p = P(P )

is equivalent to the existence of P solution of the Richards problem. The transform
P of Eq. (2.17) is

∂tp−
1

µ
∇ ·
( κ(P−1(p))

(θ′ + αP θ)(P−1(p))
K0∇p

)
− ρ0g

µ
∇ ·
(
κ(P−1(p))K0 ~e3

)
= 0.

To simplify the presentation, we introduce the notation

τ(p) =
1

µ

κ(P−1(p))

(θ′ + αP θ)(P−1(p))
.

Note that, from hypotheses (3.4)-(3.5), there exist two positive numbers τ− and τ+
such that

0 < τ− :=
κ−

µ (αP θ+ + θ′+)
≤ τ(p) ≤ τ+ :=

κ+

µαP θ−
. (3.6)

Let l = (hsoil − hbot) be the function (space depending) denoting to the total
thickness of the subsoil. We introduce the function Tl defined by

Tl(u) =
√
u+ hbot ∀u ∈ [δ2, l2],

which is extended continuously and constantly outside [δ2, l2].

Remark 3.2. To extend the solution p outside the time dependent domain Ωt, it is
necessary to impose on the function h to be less than or equal to a quantity strictly
greater than hbot. This is the reason why the small parameter δ is introduced.

Setting u = (h− hbot)
2, Equation (2.24) becomes

S0

2
∂tu−

K̃

2
∇′ · (∇′u) = −

∫ hsoil(x)

Tl(u(t,x))

∂p

∂t
dz. (3.7)

Definition. The definition of the depth h is derived from the construction of u.
Namely, for u given by (3.7), we set

h(t, x) := Tl(u). (3.8)

Remark 3.3. This definition of h allows to define the integration domain Ωt (and
thus the interface Γt) in the system (3.11)-(3.13). We emphasize that by definition,
h always remains in the interval [hbot + δ, hsoil].

We are led to consider the new problem in (u, p) completed by the boundary and
initial conditions:

S0

2
∂tu−

K̃

2
∇′ · (∇′u) = −

∫ hsoil(x)

Tl(u(t,x))

∂p

∂t
dz in (0, T )× Ωx, (3.9)

∇u · ~ν = 0 on (0, T )× ∂Ωx, u(0, x) = (h0(x)− hbot(x))2 in Ωx, (3.10)

∂tp−∇ ·
(
τ(p)K0∇p

)
− ρ0g

µ
∇ ·
(
κ(P−1(p))K0 ~e3

)
= 0 in OT , (3.11)

p
∣∣
Γt

= P(Ps) in (0, T ), ∇
(
P−1(p) + ρ0gz

)
· ~ν = 0 on (0, T )× (Γsoil ∪ Γver),

(3.12)

p(0, x, z) = P(P0)(x, z) in Ω0. (3.13)



EJDE-2022/06 DUPUIT-RICHARDS MODEL 13

Remark 3.4. Let p ∈ W (0, T ; Ω) such that p = 0 in OcT . We need to precise the

meaning of the term
∫ hsoil(x)

h(t,x)
∂p
∂t dz (h = Tl(u(t, x))):∫ hsoil(x)

h(t,x)

∂p

∂t
dz =

∫ hsoil(x)

hbot(x)

χz≥h(t,x)
∂p

∂t
dz

is the function of (H1(Ωx))′ such that ∀v ∈ H1(Ωx) ⊂ H1(Ω), for η0 > 0 small
enough 〈∫ hsoil

h(t,x)

∂p

∂t
dz, v

〉
H1(Ωx)′,H1(Ωx)

=
〈∫ hsoil

hbot

ρη0 ∗ χ{z≥(hbot+δ/2)}
∂p

∂t
dz, v

〉
H1(Ωx)′,H1(Ωx)

=
〈∂p
∂t
, ρη0

∗ χ{z≥(hbot+δ/2)} v︸ ︷︷ ︸
∈V (Ω)

〉
V ′(Ω),V (Ω)

,

where ρ ∈ C∞(R), ρ ≥ 0, with support in the unit ball such that
∫
R ρ(z)dz = 1.

We set ρη0
(z) = ρ(z/η0)/η0, η0 is chosen such that supp(ρη0

∗ χ{z≥(hbot+δ/2)}) ⊂
{z, z ≥ (hbot + δ/4)} and ρη0 ∗ χ{z≥(hbot+δ/2)} = 1 if z ≥ (hbot + 3δ/4).

The boundary condition satisfied by the unknown p at the interface Γt is then
reduced to a homogeneous Dirichlet boundary condition. We set p̄ = p − P(Ps).
Since P(Ps) is a constant, the previous system becomes

∂tp̄−∇ ·
(
τ̄(p̄))K0∇p̄

)
− ρ0g

µ
∇ ·
(
κ̄(p̄)K0 ~e3

)
= Q in OT ,

p̄
∣∣
Γt

= 0 in (0, T ),

∇
(
P−1

(
p̄+ P(Ps)

)
+ ρ0gz

)
· ~ν = 0 on (0, T )× (Γsoil ∪ Γver),

p̄(0, x, z) = P(P0)(x, z)− P(Ps) in Ω0,

where τ̄(p̄) = τ
(
p̄ + P(Ps)

)
and κ̄(p̄) = κoP−1

(
p̄ + P(Ps)

)
. We remark, that just

renaming functions τ and κ, we come back to the case P(Ps) = 0 on Γt. So, from
now, we omit the subscript ”̄.” in the previous system and we consider the system
(3.11)-(3.13) with P(Ps) = 0. The method of auxiliary domains introduced in [25]
is used. To this end, the function p is extended by zero outside the variable domain
Ωt. So, we consider the following definition of weak solution associated with system
(3.9)-(3.13).
Definition. A weak solution of (3.9)-(3.13), any function (u, p) ∈ W (0, T,Ωx) ×
W0(0, T,Ω) such that for all (φ1, φ2) ∈ L2(0, T ;H1(Ωx))× L2(0, T ;V (Ω)),∫ T

0

(S0

2
〈∂tu, φ1〉+

K̃

2

∫
Ωx

∇′u · ∇′φ1dx

+ 〈∂p
∂t
, ρη0 ∗ χ{z≥(hbot+δ/2)} φ1〉V ′(Ω),V (Ω)

)
dt = 0,

u(0, x) = (h0(x)− hbot(x))2 in Ωx,

(3.14)

∫ T

0

(
〈∂tp, φ2〉+

∫
Ω

(
τ(p)K̃0∇p+

ρ0g

µ
κ(P−1(p))K̃0 ~e3

)
· ∇φ2dx

)
dt = 0,

p = 0 in OcT , p(0, x, z) = P(P0)(x, z) in Ω0,

(3.15)



14 S. AL NAZER, C. ROSIER, M. TSEGMID EJDE-2022/06

where K̃0 = K0 in OT and K̃0 = 0 in OcT .

Theorem 3.5. Assume that there exist two real numbers θ− and κ− such that

θ(x) ≥ θ− > 0 ∀x ∈ R, κ(x) ≥ κ− > 0 ∀x ∈ R+.

Then system (3.9)-(3.13) admits a weak solution (u, p) satisfying

(a) the function u ∈ L2(0, T ;H1(Ωx)) and ∂tu ∈ L2(0, T ; (H1(Ωx))′),
(b) the function p ∈ L2(0, T ;V (Ω)) and ∂tp ∈ L2(0, T ;V ′(Ω)).

We can check that the proof of Theorem 3.5 developed below can be easily
adapted to the case of a non-zero source term Q.

Proposition 3.6. Assume that there exists δ′ > 0 such that the following a poste-
riori inequality holds

hbot + δ′ ≤ h(t, x) ≤ hsoil in (0, T )× Ωx. (3.16)

Then, the model
(
M
)

admits a weak solution (P, h) such that

(a) the function P ∈ L2(0, T ;H1(Ω)) and ∂tP ∈ L2(0, T ; (H1(Ω))′);
(b) the function h ∈ L2(0, T ;H1(Ωx)), ∂th ∈ L2(0, T ; (H1(Ωx))′).

The proof of the above proposition is a direct consequence of Theorem 3.5 as
long as the inequality 3.16 is satisfied. In this case, we turn back to the original
problem by considering the inverse transform P−1. Nevertheless, as it is done in [13,
Proposition 3], one can introduce sufficiently large “pumping/supply” source terms
allowing to control the lower and upper bounds of the solution h (and therefore to
guarantee the inequality 3.16).

4. Proof of Theorem 3.5

Let us sketch the global strategy of the proof. The problem is a strongly coupled
nonlinear system, so we apply a fixed-point approach to solve it in two steps. First,
the system is decoupled and an existence and uniqueness result for each decoupled
and linearized problem is established. The decoupled problem in p is solved by
introducing a penalized problem and by passing to the limit to come back to the
initial linearized problem. Then, we establish compactness results which allow to
prove the global existence in time of the initial problem by applying the fixed point
Schauder theorem.

4.1. Fixed point step. We now construct the framework to apply the Schauder
fixed point theorem (see [16, 34]). For the fixed point strategy, we introduce two
convex subsets (W1,W2) of W (0, T,Ωx)×W0(0, T,Ω), namely

W1 :=
{
u ∈W (0, T,Ωx) :, u(0) = u0, ‖u‖L2(0,T ;H1(Ωx)) ≤ Cu
and ‖u‖L2(0,T ;(H1(Ωx))′) ≤ C ′u

}
,

and

W2 := {p ∈W0(0, T,Ω) : p(0) = p0, ‖p‖L2(0,T ;H1(Ω)) ≤ Cp
and ‖p‖L2(0,T ;(H1(Ω))′) ≤ C ′p

}
,

constants (Cp, C
′
p) and (Cu, C

′
u) being defined thereafter.



EJDE-2022/06 DUPUIT-RICHARDS MODEL 15

Let (ū, p̄) ∈ W1 × W2, we begin by considering the unique solution u of the
linearized problem

S0

2
∂tu−

K̃

2
∇′ · (∇′u) = −

∫ hsoil(x)

h̄(t,x)

∂p̄

∂t
dz in (0, T )× Ωx, (4.1)

∇u · ~ν = 0 on (0, T )× ∂Ωx, u(0, x) = (h0(x)− hbot(x))2 in Ωx, (4.2)

where h̄(t, x) := Tl(ū(t, x)).

Remark 4.1. Note that thanks to the change of variable u = (h − hbot)
2, Equa-

tion (2.24) (which is nonlinear and degenerate in space and time) has a parabolic
structure.

Lemma 4.2. For h0 ∈ L∞(Ωx) satisfying (3.3), there exists a unique weak solution
u ∈W (0, T,Ωx) of (4.1)-(4.2) such that

‖u‖L2(0,T ;H1(Ωx)) ≤ Cu and ‖u‖L2(0,T ;(H1(Ωx))′) ≤ C ′u,

where Cu and C ′u only depend on the data of the problem.

Proof. It follows from the classical textbook [21, pp. 178-179] that for every non-
negative function ū ∈ W (0, T,Ωx) (and h̄ such that h̄(t, x) := Tl(ū(t, x))) there
exists a solution u ∈W (0, T,Ωx) of the parabolic problem with smooth coefficients

S0

2
∂tu−

K̃

2
∇′ · (∇′u) = −

∫ hsoil(x)

h̄(t,x)

∂p̄

∂t
dz,

∇u · ~ν = 0 on (0, T )× ∂Ωx, u(0, x) = (h0 − hbot)
2 in Ωx.

(4.3)

Multiplying (4.3) by u and integrating by parts over Ωx, we obtain

S0

2

d

dt

∫
Ωx

|u(t, ·)|2 dx+
K̃

2

∫
Ωx

|∇u|2 dx ≤
∣∣∣〈 ∫ hsoil

h̄(t,x)

∂p̄

∂t
dz, u

〉
H1(Ωx)′,H1(Ωx)

∣∣∣. (4.4)

Furthermore, from definition of the function Tl, we have∣∣∣〈 ∫ hsoil

h̄(t,x)

∂p̄

∂t
dz, u

〉
H1(Ωx)′,H1(Ωx)

∣∣∣
≤ ‖hsoil − hbot‖1/2∞ ‖u‖H1(Ωx)‖∂tp̄‖V ′(Ω)

≤ K̃

4

(
‖u‖2L2(Ωx) + ‖∇u‖2L2(Ωx)

)
+

1

K̃
‖hsoil − hbot‖∞‖∂tp̄‖2V ′(Ω).

Applying Gronwall’s inequality in its differential form, we obtain

‖u(t, ·)‖2L2(Ωx) ≤ e
K̃ T
2S0

(
‖u0‖2L2(Ωx) +

2

S0K̃
‖hsoil − hbot‖∞‖∂tp̄‖2L2(0,T ;V ′(Ω))

)
.

Combining the previous estimates, we deduce that

‖ u‖L2(0,T ;H1(Ωx)) ≤ C(T, S0, K̃, ‖hsoil − hbot‖∞, C ′p) := Cu.

On the other hand

‖du
dt
‖L2(0,T ;H1(Ωx)′)

= sup
‖v‖L2(0,T ;H1(Ωx))≤1

∣∣ ∫ T

0

〈du
dt
, v〉(H1(Ωx))′,H1(Ωx)dt

∣∣



16 S. AL NAZER, C. ROSIER, M. TSEGMID EJDE-2022/06

≤ 2

S0

(K̃
2
‖u(t, ·)‖L2(0,T ;H1(Ωx)) + ‖hsoil − hbot‖∞‖∂tp̄‖L2(0,T ;V ′(Ω))

)
≤ 2

S0

(K̃
2
Cu + ‖hsoil − hbot‖∞C ′p

)
:= C ′u.

The uniqueness of the solution is obvious. Indeed, if u1 and u2 are two solutions of
(4.1)-(4.2), then u = u1 − u2 satisfies

S0

2
∂tu−

K̃

2
∇′ · (∇′u) = 0 in (0, T )× Ωx,

∇u · ~ν = 0 on (0, T )× ∂Ωx, u(0, x) = 0 in Ωx.

Following the previous computations, we infer from Gronwall lemma that u = 0
a.e. in (0, T )× Ωx. This completes the proof. �

The results stated in the Lemma 3.1 require regular non-cylindrical domains in
particular with sufficiently regular boundaries (of class C1 by pieces as mentioned
by Mignot). Since in our problem, we can not guarantee as such regularity at the
interface h (which is in W (0, T,Ωx)), a regularization process is used to place our
study within the framework of Mignot [25].

We regularize h by a convolution in space. Let ψ ∈ C∞(R2), ψ ≥ 0, with
support in the unit ball such that

∫
R2 ψ(x)dx = 1. For η > 0 small enough,

we set ψη(x) = ψ(x/η)/η2. We extend h by zero outside Ωx, so we have h ∈
C([0, T ];L2(R2))∩W (0, T,R2). We define h̃ by the convolution product with respect
to the space variable

h̃ = ψη ∗ h.
Its restriction to Ωx is denoted in the same way. It fulfills h̃ ∈ C∞(Ω̄x), and as
η → 0, we have

h̃→ h strongly in C([0, T ];L2(Ωx)) ∩ L2(0, T ;H1(Ωx)).

In (3.11)-(3.13), we replace h by h̃ (the substitution appears in the space integration
domain Ωt).

Let p̄ ∈ W2 and h̃ (= ψη ∗ h) ∈ C∞(Ω̄x) where h is given by Lemma 4.2.
We consider the following linearized and regularized problem in ΩT : Find pη ∈
W0(0, T,Ω) such that for all φ ∈ L2(0, T ;V (Ω)),∫ T

0

(
〈∂tpη, φ〉+

∫
Ω

(
τ(p̄)K̃0∇pη +

ρ0g

µ
κ(P−1(p̄))K̃0 ~e3

)
· ∇φdx

)
dt = 0, (4.5)

pη = 0 in OcT and pη(0, x, z) = P(P0)(x, z) in Ω0. (4.6)

Proposition 4.3. For any η > 0, there exists a unique function pη in W0(0, T,Ω)
solution of (4.5)-(4.6). It fulfills the uniform estimates

‖pη‖L2(0,T ;H1(Ω)) ≤ Cp and ‖pη‖L2(0,T ;(H1(Ω))′) ≤ C ′p, (4.7)

where Cp and C ′p only depend on the data of the original problem (3.11)-(3.13).

Let us admit for the moment this Proposition whose the proof will be given at
the end. From now, we omit the subscript η in pη (and in uη). Let (ū, p̄) ∈W1×W2,
Lemma 4.2 and Proposition 4.3 enable to define an application F such that

W (0, T,Ωx)×W0(0, T,Ω)→W (0, T,Ωx)×W0(0, T,Ω),

F(ū, p̄) = (u, p).
(4.8)



EJDE-2022/06 DUPUIT-RICHARDS MODEL 17

The end of this subsection is devoted to the proof of the existence of a fixed point
of F in some appropriate subset. We conclude the proof of Theorem 3.5 by passing
to the limit when η → 0.

Lemma 4.4. The application F satisfies:

• There exists C a nonempty, closed, convex, bounded set in W (0, T,Ωx) ×
W0(0, T,Ω) satisfying F(C) ⊂ C,
• the application F defined by (4.8) is weakly sequentially continuous in
W (0, T,Ωx)×W0(0, T,Ω),
• there exists (u, p) ∈W1 ×W2 such that F(u, p) = (u, p).

Proof. We set C = W1 ×W2, the first point of Lemma 4.4 is obvious thanks to
Lemma 4.2 and Proposition 4.3. Indeed C is clearly a nonempty (strongly) closed
convex set in W (0, T,Ωx)×W0(0, T,Ω).

Regarding the second point of Lemma 4.4, we first note that C is compact for
the weak topology. F maps W1 ×W2 into it self. Let now (vn)n≥0 = (ūn, p̄n)n≥0

be any sequence in C which is weakly convergent in W (0, T,Ωx)×W0(0, T,Ω), and
let v = (ū, p̄) be its weak limit. We aim to show (as in [26]) that

F(vn) ⇀ F(v) in W (0, T,Ωx)×W0(0, T,Ω) as n→∞.

Since F(vn) ∈ W1 ×W2 and W1 ×W2 is weakly compact, it is sufficient to show
that there exits a subsequence (v′n) of (vn) such that F(v′n) ⇀ F(v). Extracting a
subsequence if needed we may assume without loss of generality that F(vn) ⇀ w
in W (0, T,Ωx) ×W0(0, T,Ω) as n → ∞ for some w = (u, p) ∈ W1 ×W2, and we
have to show that w and F(v) agree. Set wn = F(vn) (wn = (un, pn)), it follows
from Aubin’s Lemma that

wn → w in L2((0, T )× Ωx)× L2((0, T )× Ω) and wn(t, x)→ w(t, x) a.e.,

vn → v in L2((0, T )× Ωx)× L2((0, T )× Ω) and vn(t, x)→ v(t, x) a.e.,

∂twn ⇀ ∂tw in L2(0, T ;
(
H1(Ωx)

)′
)× L2(0, T ;V ′(Ω)),

∇wn ⇀ ∇w weakly in L2((0, T )× Ωx)× L2((0, T )× Ω).

Thanks to the Lebesgue theorem (and the properties of functions τ and Tl) we
obtain that w = F(v) (since w(0, ·) = (u(0, ·), p(0, ·)) = (u0, p0) because w ∈ C)
and the proof that F|C be weakly sequentially continuous is complete.

It follows from Schauder theorem [34] that there exists (u, p) ∈ W1 ×W2 such
that F(u, p) = (u, p). The proof of Lemma 4.4 is thus achieved.

We collect the results obtained previously. We can associate with any real num-
ber η > 0 the fixed point (uη, pη) ∈W1 ×W2 of the mapping F . It is a solution of
the system

∂tpη −∇ ·
(
τ(pη)K0∇pη

)
− ρ0g

µ
∇ ·
(
κ(P−1(pη))K0 ~e3

)
= 0 in OT , (4.9)

pη|Γt = P(Ps) in OT ,∇
(
P−1(pη) + ρ0gz

)
· ~ν = 0 on (0, T )× (Γsoil ∪ Γver),

pη(0, x, z) = P(P0)(x, z) in Ω0. (4.10)

S0

2
∂tuη −

K̃

2
∇′ · (∇′uη) = −

∫ hsoil(x)

hη(t,x)

∂pη
∂t

dz in (0, T )× Ωx, (4.11)

∇uη · ~ν = 0 on (0, T )× ∂Ωx, uη(0, x) = (h0(x)− hbot(x))2 in Ωx. (4.12)



18 S. AL NAZER, C. ROSIER, M. TSEGMID EJDE-2022/06

We can obtain similar estimates for (uη, pη) than those derived in Lemma 4.2 and
Proposition 4.3. We thus assert the existence of limit functions (extracting a sub-
sequence if needed) (u, p) ∈W (0, T,Ωx)×W0(0, T,Ω) such that

(uη, pη)→ (u, p) in L2((0, T )× Ωx)× L2((0, T )× Ω),

(uη(t, x), pη(t, x))→ (u(t, x), p(t, x)) a.e. in((0, T )× Ωx)× ((0, T )× Ω),

h̃(t, x) = ψη ∗ h(t, x)→ h(t, x) a.e. in (0, T )× Ωx,

(∂tuη, ∂tpη) ⇀ (∂tu, ∂tp) in L2(0, T ;
(
H1(Ωx)

)′
)× L2(0, T ;V ′(Ω)),

(∇uη,∇pη) ⇀ (∇u,∇p) weakly in L2((0, T )× Ωx)× L2((0, T )× Ω).

Letting η → 0 in weak formulations resulting from (4.9)-(4.12), we prove the exis-
tence of a weak solution (u, p) of problem (3.9)-(3.13). This completes the proof of
Theorem 3.5.

4.2. Proof of Proposition 4.3. Again, we omit the subscript η in pη. The proof
of Proposition 4.3 is done by introducing a penalized problem and by passing to
the limit to come back to the linearized problem (4.5)-(4.6). We thus consider the
weak solution p of the linearized problem (4.5)-(4.6). So we look for p ∈W0(0, T,Ω)
such that, for all φ ∈ L2(0, T ;V (Ω)),∫ T

0

(
〈∂tp, φ〉+

∫
Ω

(
τ(p̄)K̃0∇p+

ρ0g

µ
κ(P−1(p̄))K̃0 ~e3

)
· ∇φ

)
dt = 0.

We first remark that the solution of system (4.5)-(4.6) is unique. Indeed, if p1 and
p2 are two solutions of (4.5)-(4.6), then q = p1 − p2 satisfies∫ T

0

(
〈∂tq, φ〉+

∫
Ωt

τ(p̄)K0∇q · ∇φdx
)
dt = 0.

Then, taking φ = q and using the fourth point of Lemma 3.1, we conclude that

1

2

∫
ΩT

q2(T, x) dx +

∫ T

0

∫
Ωt

τ(p̄)K0|∇q|2 dx dt = 0,

since q(0, ·) = 0. We infer from this equality that q = 0 a.e. in (0, T ) × Ω (since
q = 0 on the interface Γt).

We will define a family of approximate problems which are linear parabolic prob-
lems in the cylindrical domain (0, T )×Ω, and whose the solution restricted to the
set OT will converge to the solution p of the linearized equation (4.5).

Step 1. Penalized problems. Let ε > 0, we now consider the following penalized
problem on Ω: Find pε ∈W0(0, T,Ω) such that for all φ ∈ L2(0, T ;D(Ω̄)) null in a
neighborhood of Γbot,∫ T

0

〈∂tpε, φ〉 dt+

∫ T

0

∫
Ω

(
K̃0τ(p̄)∇pε +

ρ0g

µ
κ(P−1(p̄))K̃0 ~e3

)
· ∇φdx dt

+

∫
OcT
∇pε · ∇φdx dt+

1

ε

∫
OcT

pεφdx dt = 0,

(4.13)

pε(0, x, z) = P(P0)(x, z) in Ω0 and pε(0, x, z) = 0 in Ω \ Ω0. (4.14)

We aim to state that the penalized system (4.13)-(4.14) admits a unique solution
pε which tends to the solution of problem (4.5)-(4.6) when ε→ 0. Equation (4.13)



EJDE-2022/06 DUPUIT-RICHARDS MODEL 19

can be written, for all φ ∈W0(0, T,Ω),∫ T

0

∫
Ω

τ(p̄)K̃0∇pε · ∇φdx dt+

∫
OcT
∇pε · ∇φdx dt+

1

ε

∫
OcT

pεφdx dt︸ ︷︷ ︸
Aε(,pε,φ)

+

∫ T

0

〈∂tpε, φ〉 dt

= −
∫ T

0

∫
Ω

ρ0g

µ
κ(P−1(p̄))K̃0 ~e3 · ∇φ

)
dx dt︸ ︷︷ ︸

Lε(φ)

.

(4.15)

From (3.6), we establish that the coefficients of Aε are in L∞((0, T )×Ω). Moreover,
we have

Aε(p, p) ≥ inf(1,K0τ−,
1

ε
)‖p‖L2(0,T ;V (Ω)), ∀p ∈ L2(0, T ;V (Ω)).

We directly check that Lε is a linear form on L2(0, T ;V (Ω)). We thus deduce the
existence and uniqueness for the system (4.13)-(4.14).

Step 2. Limit when ε→ 0. We first derive some uniform estimates with respect
to ε (and η). Multiplying (4.15) by pε and integrating by parts over Ω, we obtain
for all s ≤ T that∫ s

0

〈∂tpε, pε〉 dt+

∫ s

0

∫
Ω

τ(p̄)K̃0|∇pε|2 dx dt+

∫
Ocs
|∇pε|2 dx dt+

1

ε

∫
Ocs
p2
ε dx dt︸ ︷︷ ︸

I1

= −
∫ s

0

∫
Ω

ρ0g

µ
κ(P−1(p̄))K̃0 ~e3 · ∇pε dx dt︸ ︷︷ ︸

I2

, ∀φ ∈W0(0, T,Ω).

Then, applying Lemma 3.1 to the first term, we obtain∣∣ ∫ s

0

〈∂tpε, pε〉V ′(Ω),V (Ω) dt
∣∣ =

1

2

(∫
Ω

p2
ε(s, ·)dx−

∫
Ω

p2
ε(0, ·)dx

)
.

Also

|I1| ≥ K0 τ−‖∇pε‖2L2(Os) + ‖∇pε‖2L2(Ocs) +
1

ε

∫
Ocs
p2
ε dx dt,

|I2| ≤
ε1
2
‖∇pε‖2L2(Os) +

1

2 ε1
T meas(Ω)

(ρ0gκ+K0

µ

)2

.

By taking ε1 = K0 τ−, we deduce directly from these estimates that the sequence
{pε} is bounded in L2(0, T ;V (Ω)) and the sequence {pε/

√
ε} is bounded in L2(OcT ).

Indeed, applying Gronwall’s Lemma, there exists a constant Cp depending only on
the data such that

‖pε‖2L2(0,T ;V (Ω)) ≤ Cp. (4.16)

More precisely, we have

min
(
1,
K0τ−

2

)
‖∇pε‖2L2((0,T )×Ω) ≤

T

2K0 τ−
meas(Ω)

(ρ0gκ+K0

µ

)2

+
1

2

∫
Ω

p2
0dx.



20 S. AL NAZER, C. ROSIER, M. TSEGMID EJDE-2022/06

We thus can extract subsequences {pε}, { 1√
ε
pε} (not relabeled for convenience) and

there exist q ∈ L2(0, T ;V (Ω)) and q′ ∈ L2(OcT ) such that

pε ⇀ q weakly in L2((0, T )× Ω) (4.17)

1√
ε
pε ⇀ q′ weakly in L2(OcT ) (4.18)

∇pε ⇀ ∇p weakly in L2((0, T )× Ω). (4.19)

It results from the first two convergences that

pε|O′T ⇀ q|O′T weakly in L2(OcT )

pε|O′T =
√
ε× 1√

ε
pε ⇀ 0 weakly in L2(OcT ),

so

q|OcT = 0. (4.20)

Moreover, since q ∈ L2(0, T ;V (Ω)), we infer from the second result of Lemma
3.1, that we can define γ(q) on Γt. Furthermore thanks to (4.20), γ(q) = 0 on
Γt, 0 ≤ t ≤ T , and thus q|OT ∈ F (OT ). We must now check that q|OT satisfies
(4.5).

It remains to be established that

Dt(q|OT ) ∈ F ′(OT ) and Dt(pε|OT ) ⇀ Dt(q|OT ) in F ′(OT ).

Let ψ ∈ BT (OT ), the function ψ is extended by zero on OcT , the extension denoted

ψ̃ thus belongs to W0(0, T,Ω). Hence taking φ = ψ̃ in (4.13) leads to∫ T

0

〈∂tpε, ψ̃〉 dt+

∫ T

0

∫
Ω

(
K̃0τ(p̄)∇pε +

ρ0g

µ
κ(P−1(p̄))K̃0 ~e3

)
· ∇ψ̃ dx dt = 0,

that we can write as follows

−
∫ T

0

〈∂tψ̃, pε〉 dt+

∫ T

0

∫
Ω

(
K̃0τ(p̄)∇pε +

ρ0g

µ
κ(P−1(p̄))K̃0 ~e3

)
· ∇ψ̃ dx dt

=

∫
Ω0

p0(x)ψ0(0,x) dx.

By letting ε→ 0, we obtain

−
∫ T

0

〈∂tψ, q|OT 〉 dt+

∫ T

0

∫
Ωt

(
K0τ(p̄)∇q|OT +

ρ0g

µ
κ(P−1(p̄))K0 ~e3

)
· ∇ψ dx dt

=

∫
Ω0

p0(x)ψ0(0,x) dx.

From this equation, we deduce that, taking ψ ∈ D(ŌT ) null in a neighborhood of
Γt

Dt(q|OT )−∇ ·
(
K0τ(p̄)∇q|OT +

ρ0g

µ
κ(P−1(p̄))K0 ~e3

)
= 0.

Since F (OT ) is the closure inH0,1(OT ) of functions ofD(ŌT ) null in a neighborhood
of Γt, the previous equality holds true in F ′(OT ). Moreover the solution pε|OT of
(4.15) verifies

Dt(pε|OT )−∇ ·
(
K0τ(p̄)∇pε|OT +

ρ0g

µ
κ(P−1(p̄))K0 ~e3

)
= 0 in F ′(OT ).



EJDE-2022/06 DUPUIT-RICHARDS MODEL 21

Letting ε→ 0 and we infer from convergences (4.17)-(4.19) that

∇ ·
(
τ(p̄)∇pε|OT +

ρ0g

µ
κ(P−1(p̄))~e3

)
︸ ︷︷ ︸

=Dt(pε|OT )

⇀ ∇ ·
(
τ(p̄)∇q|OT +

ρ0g

µ
κ(P−1(p̄))~e3

)
︸ ︷︷ ︸

=Dt(q|OT )

,

weakly in F ′(OT ). Thus q|OT is the unique solution of (4.5)-(4.6), and the limit of
pε|OT being independent of the chosen subsequence, the whole sequence converges
towards q|OT . Moreover, the first inequality in (4.7) is obtained for the solution
q ∈ L2(0, T ;V (Ω)) of system (4.5)-(4.6) in the same way as for the estimate (4.16)
obtained for pε. Finally, as was done in Lemma 4.2, we deduce from the first
inequality of (4.7) that

‖∂tq‖2L2(0,T,V ′(Ω)) ≤ C
′
p,

where C ′p depends on the data and on Cp. This completes the proof of Proposition
4.3.

References

[1] M. B. Abbott, J. C. Bathurst, J. A. Cunge, P. E. O’connell, J. Rasmussen; An introduction

to the european hydrological system - systeme hydrologique europeen, “she”, 2: Structure of
a physically-based, distributed modelling system. Journal of Hydrology, 87(1) (1986), 61–77.

[2] P. Ackerer, A. Younes; Efficient approximations for the the simulation of density driven flow
in porous media, Adv. Water Res., Vol. 31 (2008), 15–27.

[3] H. W. Alt, S. Luckhaus; Quasilinear elliptic-parabolic differential equations, Math. Z., Vol.

1 (1983), 311–341.
[4] H. W. Alt, E. Di Benedetto; Nonsteady flow of water and oil through inhomogeneous porous

media, Ann. Scuola Norm. Sup. Pisa , 12 (1985), 335–392.

[5] J. Bear; Dynamics of fluids in porous media, Elsevier, New-York (1972).
[6] J. Bear, A. Verruijt; Modeling groundwater flow and pollution, Springer, Netherlands (1987).

[7] H. Berninger, M. Ohlberger, O. Sander, K. Smetana; Unsaturated subsurface flow with sur-

face water and nonlinear in- and outflow conditions. Mathematical Models and Methods in
Applied Sciences, 24(05) (2014), 901–936.

[8] C. Bourel, C. Choquet, C. Rosier, M. Tsegmid; Modelling of shallow aquifers in interaction

with overland water, Appl. Math. Model., 81 (2020), 727–751.
[9] R. H. Brooks, A. T. Corey; Hydraulic Properties of Porous Media, Colorado State University

Hydrology Papers, (1964).
[10] X. Chen, A. Friedman, T Kimura; Nonstationary filtration in partially saturated porous

medium, Euro. J. Appl. , 5 (1994), 405–429.
[11] C. Choquet,, M. M. Diédhiou, C. Rosier; Derivation of a Sharp-Diffuse Interfaces Model for

Seawater Intrusion in a Free Aquifer. Numerical Simulations, SIAM J. Appl. Math., 76
(2016), no. 1, 138–158.

[12] C. Choquet, M. M. Diédhiou, C. Rosier; Mathematical analysis of a sharp-diffuse interfaces
model for seawater intrusion, J. Differential Equations, 259 (2015), no. 8, 3803–3824.

[13] C. Choquet, C. Rosier, L. Rosier; Well posedness of general cross-diffusion systems, J. Dif-
ferential Equations, 300 (2021), 386–425.

[14] H. Darcy; Les fontaines publiques de la ville de Dijon; exposition et application des principes
à employer dans les questions de distribution d’eau, Victor Dalmont, Editeur, Paris, (1856).

[15] J. Dupuit; Etudes thd́oriques et pratiques sur le mouvement des eaux dans les canaux
dćouverts et à travers les terrains permśables, 2ème édition, Dunod, Paris, (1863).

[16] L. C. Evans; Partial differential equations, American Mathematical Society (1998).

[17] C. W. Fetter; Hydrogeology: A shot history, Part 2, Ground Water, 42 (2004), 949–953.
[18] J. Hulshof, N. Wolanski; Monotone flows in N-dimensional partially saturated porous media:

Lipschitz continuity of the interface, Arch. Rat. Mech. Anal., 102 (1988), 287–305.

[19] H. Kong, R.Lian; Free boundary value problem for compressible magnetohydrodynamic equa-
tions, Electron. J. Differential Equations, 2020 (2020), no. 11, 1–16.



22 S. AL NAZER, C. ROSIER, M. TSEGMID EJDE-2022/06

[20] J. Kong, P. Xin, Z. Yao Song, L. Li; A new model for coupling surface and subsurface water

flows: With an application to a lagoon, Journal of Hydrology, 390(1) (2010), 116–120.

[21] O. A. Ladyzhenskaja, V. A. Solonnikov, N. N. Ural’ceva; Linear and quasilinear equations
of parabolic type, Amer. Math. Soc. (1968).

[22] J. L. Lions; Sur les problèmes mixtes pour certains systèmes paraboliques dans des ouverts

non cylindriques, Annales de l’Institut Fourier, 7 (1957), 143–182.
[23] J. L. Lions, E. Magenes; Problèmes aux limites non homogènes, Vol. 1, Dunod (1968).

[24] C. M. Marle, Henry Darcy et les ’ecoulements de fluides en milieux poreux, Oil and Gas

science and technology Re. IFP, Vol. 61 (5) (2006), 599–609.
[25] A. L. Mignot; Méthodes d’approximation des solutions de certains problèmes aux limites

linéaires; Rendiconti del Seminario Matematico della Università di Padova, tome 40 (1968),

p. 1-138.
[26] K. Najib, C. Rosier; On the global existence for a degenerate elliptic-parabolic seawater

intrusion problem, Math. Comput. Simulation, 81 Issue 1 (2011), 2282–2295.
[27] M. F. Pikul, R. L. Street, I. Remson; A numerical model based on coupled one-dimensional

richards and boussinesq equations, Water Resources Research, 10(2) (1974), 295–302.

[28] R. E. Showalter, N. Su; Partially saturated flow in a poroelastic medium, Disc. Cont. Dyn.
Syst. Ser. B (2001), 403-420.

[29] J. Simon; Compact sets in the space Lp(0, T, B), Ann. Mat. Pura Appl., 146 (4) (1987),

65–96.
[30] P. Sochala, A. Ern, S. Piperno; Mass conservative bdf-discontinuous galerkin/explicit finite

volume schemes for coupling subsurface and overland flows, Computer Methods in Applied

Mechanics and Engineering, 198(27) ( 2009), 2122–2136.
[31] M. Tsegmid; Modélisation d’aquifères peu profonds en interaction avec les eaux de surfaces,

Phd Thesis, 2019 ULCO. www.theses.fr/2019DUNK0526

[32] M. T. Van Genuchten; A closed-form equation for predicting the hydraulic conductivity of
unsaturated soils 1, Soil science society of America journal, 44(5) (1980), 892–898.

[33] H. M. Yin, A singular-degenerate free boundary problem arising from the moisture evapora-
tion in a partially saturated porous medium, Ann. Mat. Pura Appl. , 161 (1992), 379-397.

[34] E. Zeidler, Nonlinear functional analysis and its applications, Part 1, Springer Verlag, (1986).

Safaa Al Nazer

Université du Littoral Côte d’Opale, UR 2597, LMPA, Laboratoire de Mathématiques
Pures et Appliquées Joseph Liouville, F-62100 Calais, France

Email address: safaaalnazer9@gmail.com

Carole Rosier

Université du Littoral Côte d’Opale, UR 2597, LMPA, Laboratoire de Mathématiques

Pures et Appliquées Joseph Liouville, F-62100 Calais, France
Email address: Carole.Rosier@univ-littoral.fr

Munkhgerel Tsegmid
School of Applied Sciences, Mongolian University of Science and Technology, Ulaan-
baatar, Mongolia

Email address: tsemuugii@must.edu.mn


	1. Introduction
	2. Derivation of the model
	2.1. Geometry
	2.2. Conservation laws
	2.3. State equation for the fluid compressibility
	2.4. Permeability tensor K0
	2.5. Hypothesis
	2.6. Model coupling vertical 3d-Richards flow and Dupuit horizontal flow

	3. Mathematical setting and main results
	3.1. Notation and auxiliary results
	3.2. Main results

	4. Proof of Theorem 3.5
	4.1. Fixed point step
	4.2. Proof of Proposition 4.3

	References