MATHEMATICAL ANALYSIS OF A DUPUIT-RICHARDS MODEL

. This article is concerned with a model which is an alternative to the 3D-Richards equation to describe the ﬂow of water in shallow aquifers. The model couples the two dominant types of ﬂow existing in the aquifer. The ﬁrst is described by the classic Richards problem in the upper capillary fringe. The second results from Dupuit’s approximation after vertical integration of the conservation laws between the bottom of the aquifer and the saturation interface. The ﬁnal model consists of a strongly coupled system of parabolic-type pde which are deﬁned in a time-dependent domain. First, we show how taking the low compressibility of the ﬂuid 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.


Introduction
Populated areas are increasingly affected by contamination of soil and groundwater.Many modeling approaches are 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 involved processes (chemical, hydrogeological, anthropic, ...) acting in a wide range of temporal and geometrical length scales.But we can 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.Indeed, 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).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 paper, a model belonging to the "Dupuit-Richards" model class is analyzed.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 conductivity tends towards zero in the unsaturated part of the aquifer.Finally, the coupling of the 1D Richards equation with the 2D Dupuit approximation numerically 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 nonlinearities, 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 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 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 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 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 nonlinearity 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 quasilinear 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 conservation 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.
The document 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.

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 = (x 1 , x 2 ) ∈ R 2 , z ∈ R, the usual coordinates.
2.1.Geometry.The aquifer is represented by a three-dimensional domain Ω := Ω x × (h bot , h soil ), Ω x ⊂ R n with n ≥ 2 (x = (x 1 , x 2 )), function h bot (respect.h soil ) describing its lower (respect.upper) topography.The upper and lower surfaces are thus defined by the graph of the functions h bot = h bot (x) and h soil = h soil (x), x ∈ Ω x .We assume that More precisely the domain is given by: We always denote by ν the outward unit normal and e 3 is the unitary vertical vector pointing up.The boundary ∂Ω of Ω is divided into three zones (bottom, top and vertical) The model divides the flow description into two subregions of Ω (possibly timedependent) 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 h bot ≤ h ≤ h soil , we introduce: and Ω t := (x, z) ∈ Ω | z > h(t, x) , (2.3) (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 where ρ and µ are respectively the density and the viscosity of the fluid, K 0 is the permeability of the soil, κ(P ) is the relative conductivity and g the gravitational acceleration constant.Introducing the hydraulic head H defined by the previous equation is written as follows: 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 following equation where Q denotes a generic source term (for production and replenishment).The function θ is the volumetric moisture content defined by 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 latter.

2.3.
State equation for the fluid compressibility.The fluid is considered compressible by assuming that pressure P is related to the density ρ as follows (cf.[11]): The real number α P ≥ 0 represents the fluid compressibility coefficient and P 0 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 K 0 .The nonlinear hydraulic conductivity K is given by K = κ(P ) ρ0g µ K 0 .The soil transmission properties are characterized by the porosity function φ and the permeability tensor K 0 (x, z).The matrix K 0 is a 3 × 3 symmetric positive definite tensor which describes the conductivity of the saturated soil at the position (x, z) ∈ Ω.We introduce (2.9) 2.5.Hypothesis.Let us now list the assumptions on the fluid and medium characteristics 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 get: Second consequence is ∇ρ • q 1 which leads to the following simplification in the mass conservation equation (2.7): 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 ∂ t P in the latter equation, we get After simplification by ρ > 0, we finally obtain (2.12) Equivalently, using the hydraulic head (2.5) and the Darcy law (2.11),(2.12) can be written We notice that if the fluid is assumed incompressible, α P = 0, then Eq. (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 P s which is a fixed real number.The fully-saturated part of the medium corresponds to the region {x; P (., x) > P s }, while it is partially-saturated in the capillary fringe {x; P d < P (., x) ≤ P s }.The dry part is defined by the set {x; P (., x) ≤ P d }.The moisture content is such that ) where s 0 > 0 corresponds to a residual saturation which is positive.The associated relative hydraulic mobility is then defined by (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 dependance 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 ≥ P s κ(P ) = 1 ⇐⇒ P ≥ P s . (2.16) In particular, the water pressure is greater than the bubbling pressure P s 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 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 (2.17) The effective velocity q is given by We emphasize that the model (2.17) depends by definition on the depth h which is expected to belong to the interval (h bot , h soil ).
• 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 h bot and h.Since θ(P ) = φ in the saturated zone, the vertical average (2.13) leads to denotes the thickness of the saturated zone and Q the source term representing distributed surface supply of fresh water into the free aquifer: Applying Leibnitz rule to the first term in the left-hand side yields: We denote by H the vertically averaged hydraulic head (we remind that κ(P ) = 1 for z ∈ (h bot , h)).The averaged mass conservation law for the freshwater in the saturated zone finally reads In this equation, the term B f 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 (x 1 , x 2 , z, t) = 0 ⇔ z − h(x 1 , x 2 , 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 Approximation of the flux q |z=h + • ∇(z − h): The flux q |z=h + • ∇(z − h) expresses mass transfers between the two parts of the aquifer.As it is done in [8], we approximate the flux by This approximation comes from the hypothesis of an almost null horizontal hydraulic conductivity (i.e.K xx (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 h soil 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 = h soil : Since the lower layer is impermeable, there is no flux across the boundary z = h bot : Besides, the pressure is continuous across Γ t , it follows that with K 0 ρ 0 g µ dz and S 0 = ρ 0 g φ α P . (2.26) The homogeneous Neumann condition on ∂Ω x is assumed to simplify the presentation.
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 following 3d-Richards equation holds where the saturation pressure P s is assumed to be constant with respect to the time and to the space.The effective velocity q is given by • In Ω − t the pressure P satisfies in Ω x .

Mathematical setting and main results
The problem (M) being a problem with free boundary, we are going to define the general framework of parabolic equation in non cylindrical domain, 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.
where h is the position of the interface Γ t .We set We define It is an Hilbert space endowed with the norm F (O T ) denotes the closure in H 0,1 (O T ) of functions of D( ŌT ) null in a neighborhood of Γ t and F (O T ) its topological dual.Besides, we introduce endowed with the Hilbertian norms We now state some auxiliary results proved in [22] Lemma 3.1.If O T is sufficiently regular, we thus have where For the sake of brevity, we shall write H 1 (Ω) = W 1,2 (Ω) and are dense and compact.For any T > 0, let W 0 (0, T, Ω) denotes the space . The following embeddings are continuous ([23] prop.2.1 and thm 3.1, chapter 1) is compact (Aubin's Lemma, see [29]).
In the same way, we introduce W (0, T, Ω x ) the space endowed with the Hilbertian norm . The same compacity results hold true in this case too.

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 K 0 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 P 0 ∈ H 2 (Ω) satisfies the compatibility condition P 0 (x, h 0 ) = P s in Ω 0 .Let δ ∈ R be a positive number, we assume that h 0 ∈ L ∞ (Ω x ) is such that Functions θ and κ are pressure-dependent and we assume (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 Eq. (2.17A direct computation gives P (P ) = θ (P ) + α P θ(P ) > α P θ − > 0, indeed by previous hypothesis, we have θ (P ) ≥ 0 and θ(P ) > φs 0 .Since θ ∈ C 1 (R), there exists θ + > 0 such that 0 ≤ θ ≤ θ + on the interval [P d , P s ].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 To simplify the presentation, we introduce the notation .
Note that, due to hypotheses (3.4)-(3.5),there exist two positive reals τ − and τ + such that 0 Let l = (h soil − h bot ) be the function (space depending) denoting to the total thickness of the subsoil.We introduce the function T l defined by which is extended continuously and constantly outside [δ 2 , l 2 ].
Remark 3.2.In order 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 h bot .This is the reason why the small parameter δ is introduced.
Setting u = (h − h bot ) 2 , Eq. (2.24) becomes 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) := T l (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.12).We emphasize that by definition, h always remains in the interval We are led to consider the new problem in (u, p) completed by the boundary and initial conditions: p Γt = P(P s ) in (0, T ), ∇ P −1 (p) + ρ 0 gz • ν = 0 on (0, T ) × (Γ soil ∪ Γ ver ), p(0, x, z) = P(P 0 )(x, z) in Ω 0 . (3.12) Remark 3.4.Let p ∈ W (0, T ; Ω) such that p = 0 in O c T .We need to precise the meaning of the term 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(P s ).Since P(P s ) is a constant, the previous system becomes p Γt = 0 in (0, T ), ∇ P −1 p + P(P s ) + ρ 0 gz • ν = 0 on (0, T ) × (Γ soil ∪ Γ ver ), p(0, x, z) = P(P 0 )(x, z) − P(P s ) in Ω 0 , where τ (p) = τ p + P(P s ) and κ(p) = κoP −1 p + P(P s ) .We thus remark, that just renaming functions τ and κ, we come back to the case P(P s ) = 0 on Γ t .So, from now, we omit the subscript "." in the previous system and we consider the system (3.11)-(3.12)with P(P s ) = 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 the system (3.9)-(3.12): Definition: We call weak solution of problem (3.9)-(3.12),any solution where K0 = K 0 in O T and K0 = 0 in O c T .We state and prove the following existence result: Theorem 3.5.Assume that there exist two real numbers θ − and κ − such that The system (3.9)-(3.12)admits a weak solution (u, p) satisfying (a) the function u ∈ L 2 (0, T ; We can check that the proof of Theorem 3.5 developed below may 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 posteriori inequality holds Then, the model M admits a weak solution (P, h) such that (a) the function P ∈ L 2 (0, T ; H 1 (Ω)) and The proof of Proposition 3.6 is a direct consequence of Theorem 3.5 as soon as the inequality 3.15 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 Proposition 3 of [13], 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.15).Next section is devoted to the proof of Theorem 3.5.

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 (W 1 , W 2 ) of W (0, T, Ω x ) × W 0 (0, T, Ω), namely constants (C p , C p ) and (C u , C u ) being defined thereafter.Let (ū, p) ∈ W 1 ×W 2 , we begin by considering the unique solution u of the following linearized problem where h(t, x) := T l (ū(t, x)).
Remark 4.1.Note that thanks to the change of variable u = (h − h bot ) 2 , Eq.
(2.24) (which is nonlinear and degenerate in space and time ) has now a parabolic structure.
3), there exists a unique weak solution u ∈ W (0, T, Ω x ) of (4.1)-( 4.2) such that where C u 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 nonnegative function ū ∈ W (0, T, Ω x ) (and h s.t.h(t, x) := T l (ū(t, x))) there exists a solution u ∈ W (0, T, Ω x ) of the parabolic problem with smooth coefficients Multiplying Eq. (4.3) by u and integrating by parts over Ω x , we obtain Furthermore, by definition of the function T l , we have Applying Gronwall's inequality in its differential form, we get Combining the previous estimates, we deduce that On the other hand The uniqueness of the solution is obvious.Indeed, if u 1 and u 2 are two solutions of (4.1)-(4.2),then u = u 1 − u 2 satisfies Following the previous computations, we infer from Gronwall lemma that u = 0 a.e. in (0, T ) × Ω x .This ends the proof of Lemma 4.2.
The results stated in the lemma 3.1 require regular non-cylindrical domains in particular with sufficiently regular boundaries (of class C 1 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 thus regularize h by convolution in space.Let ψ ∈ C ∞ (R 2 ), ψ ≥ 0, with support in the unit ball such that R 2 ψ(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 ]; L 2 (R 2 ))∩W (0, T, R 2 ).We define h by the convolution product with respect to the space variable h = ψ η * h.
We thus consider the following linearized and regularized problem in Ω T : Find 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) ∈ W 1 ×W 2 , Lemma 4.2 and Proposition 4.3 enable to define an application F such that: The end of the present 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.
Proof.We set C = W 1 × W 2 , 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 ) × W 0 (0, T, Ω).Regarding the second point of Lemma 4.4, we first note that C is compact for the weak topology.F maps W 1 × W 2 into it self.Let now (v n ) n≥0 = (ū n , pn ) n≥0 be any sequence in C which is weakly convergent in W (0, T, Ω x ) × W 0 (0, T, Ω), and let v = (ū, p) be its weak limit.We aim to show (as in [26]) that Extracting a subsequence if needed we may assume without loss of generality that F(v n ) w in W (0, T, Ω x ) × W 0 (0, T, Ω) as n → ∞ for some w = (u, p) ∈ W 1 × W 2 , and we have to show that w and F(v) agree.Set w n = F(v n ) (w n = (u n , p n )), it follows from Aubin's Lemma that Thanks to Lebesgue theorem (and the properties of functions τ and T l ) we obtain that w = F(v) (since w(0, .)= (u(0, .),p(0, .))= (u 0 , p 0 ) 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) ∈ W 1 × W 2 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 number η > 0 the fixed point (u η , p η ) ∈ W 1 × W 2 of the mapping F. It is a solution of the system : ) 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 subsequence if needed) (u, p) ∈ W (0, T, Ω x ) × W 0 (0, T, Ω) such that Letting η → 0 in weak formulations resulting from (4. We first remark that the solution of system (4.4)-(4.5) is unique.Indeed, if p 1 and p 2 are two solutions of (4.4)-(4.5),then q = p 1 − p 2 satisfies Then, taking φ = q and using the fourth point of Lemma 3.1, we conclude that 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 problems in the cylindrical domain (0, T ) × Ω, and whose the solution restricted to the set O T will converge to the solution p of the linearized equation (4.4).
Step 1. Penalized problems Let > 0, we now consider the following penalized problem on Ω: Find p ∈ W 0 (0, T, Ω) s.We aim to state that the penalized system (4.12)-(4.13)admits a unique solution p which tends to the solution of problem (4.4)-(4.5)when → 0. Eq. (4.12) can be written, ∀φ ∈ W 0 (0, T, Ω) We thus can extract subsequences {p }, { 1 √ p } (not relabeled for convenience) and there exist q ∈ L 2 (0, T ; V (Ω)) and q ∈ L 2 (O c T ) such that p q weakly in L 2 ((0, T ) × Ω) (4.16) 1 √ p q weakly in L 2 (O c T ) (4.17Moreover, since q ∈ L 2 (0, T ; V (Ω)), we infer from the second result of Lemma 3.1, that we can define γ(q) on Γ t .Furthermore thanks to (4.19), γ(q) = 0 on Γ t , 0 ≤ t ≤ T , and thus q| O T ∈ F (O T ).We must now check that q| O T satisfies (4., weakly in F (O T ).Thus q|O T is the unique solution of (4.4)-(4.5),and the limit of p | O T being independent of the chosen subsequence, the whole sequence converges towards q|O T .Moreover, the first inequality in (4.6) is obtained for the solution q ∈ L 2 (0, T ; V (Ω)) of system (4.4)-(4.5) in the same way as for the estimate (4.15) obtained for p .Finally, as was done in Lemma 4.2, we deduce from the first inequality of (4.6) that where C p depends on the data and on C p .This ends the proof of Proposition 4.3.

3. 1 .
Notations and auxiliary results.For any T > 0, let O T be the open domain of R + × Ω defined by ), namely assumptions (3.4)-(3.5)are sufficient to define the primitive function P such that P(P ) = θ(P ) + α P P θ(s) ds.