A degenerate elliptic-parabolic system arising in competitive contaminant transport
Margarida Baía, Farid Bozorgnia, Léonard Monsaingeon, Juha Videman
AA DEGENERATE ELLIPTIC-PARABOLIC SYSTEM ARISING INCOMPETITIVE CONTAMINANT TRANSPORT
MARGARIDA BAÍA, FARID BOZORGNIA, LÉONARD MONSAINGEON, AND JUHA VIDEMAN
Abstract.
In this work we investigate a coupled system of degenerate and non-linear partial differential equations governing the transport of reactive solutes ingroundwater. We show that the system admits a unique weak solution provided thenonlinear adsorption isotherm associated with the reaction process satisfies certainphysically reasonable structural conditions. We conclude, moreover, that the soluteconcentrations stay non-negative if the source term is componentwise non-negativeand investigate numerically the finite speed of propagation of compactly supportedinitial concentrations, in a two-component test case. Introduction
Contaminant transport in groundwater and subsurface environments is a complexprocess given the combination of phenomena, phases and interfaces present in suchenvironments. The transport can be influenced by an array of mechanical and chemicalinteractions, between or among the different constituent phases and species, such asadvection, diffusion, dispersion and sorption, cf. [6, 7].In this work, we will study a system of partial differential equations modeling multi-component contaminant transport (transport of several pollutants). For expediency,we focus on the case where the transport is dominated by adsorption and diffusion.The adsorption process, adherence of a pollutant on the solid matrix at the fluid-solidinterface, is assumed to be in equilibrium and of Freundlich type. In Section 2, we showthat these assumptions lead to the initial-boundary value problem(1.1) ∂ t b ( u ) − ∆ u = f in (0 , T ) × Ω ,u = u b on (0 , T ) × Γ ,b ( u | t =0 ) = b in Ω , for the unknown solute concentrations u : [0 , T ] × Ω → R m , where f : [0 , T ] × Ω → R m , b : Ω → R m and u b : [0 , T ] × Γ → R m are given, Ω ⊂ R d is a smooth domain where thetransport process takes place, Γ = ∂ Ω , T > and m ≥ . Moreover, the vector field b : R m → R m , associated with the Freundlich adsorption isotherm, can be written as agradient of a non-negative convex function and is such that D u b ( u ) → ∞ when u → ,i.e., the elliptic-parabolic system (1.1) becomes singular at u = 0 . The system may thusexhibit finite speed of propagation of compactly supported initial concentrations givingrise to free boundaries that separate the regions where the solute concentrations vanishfrom those where they are positive.Problem (1.1) falls under the general quasi-linear elliptic-parabolic systems addressedby Alt and Luckhaus in their seminal paper [2]. Even if our method of proof is largelyinfluenced by these results, we wish to present here a detailed, and at the same timemore straightforward, existence proof. Our approach, based on Rothe’s method ([18])and on solving a convex minimization problem at each time step, provides also a simplemethod for the numerical approximation of solutions to system (1.1).If the nonlinear term b behind the time derivative in (1.1) is invertible, as a func-tion from R m to R m , one can make a change of variables and pass the degeneracy tothe diffusion term. In that case, system (1.1) reduces to a coupled system of nonlineardiffusion equations of porous medium type or, in the one-component case, to a general-ized porous medium equation. This inversion is possible for most adsorption isotherms,thus in the one-component case the existence theory follows from nowadays well-knownresults on nonlinear diffusion equations (cf. [24, 8]). We note, however, that most nu-merical methods applied to the single species case keep the equation in the form (1.1),see [10, 9, 4, 5, 1]. a r X i v : . [ m a t h . A P ] J a n BAÍA, BOZORGNIA, MONSAINGEON, AND VIDEMAN
The multi-component case is much trickier, as systems of nonlinear PDEs often are,and the only existence proof known to us leans heavily on a special structure of thesystem wherein the sum of the solute concentrations solve a scalar generalized porousmedium equation, see [14]. Our system (1.1) does not possess such structure but onthe other hand b satisfies the assumption of being a gradient of scalar function whichthe system studied in [14] did not, so the existence results are mutually exclusive. Atthe same time, we note that our numerical results for the two-component case arequalitatively very similar to the ones presented in [14].The paper is organized as follows. In Section 2, we present an overview of the physicalmodeling of competitive contaminant transport in groundwater and derive our modelproblem. In Section 3, we introduce our notation and some auxiliary results needed forthe proof of our main existence result that is proved in Section 4, where we start by givinga mathematically precise meaning to the solutions we are looking for and we introducethe structural assumptions imposed on the nonlinear isotherm b . In Section 5, making anadditional assumption on the isotherm potential φ , we show that the solution stays non-negative provided the source data is componentwise non-negative which of course is to beexpected on physical grounds. Finally, in Section 6, we present a simple finite-differencescheme which mimics our existence proof. We solve a two-component test problem andobserve that the scheme captures quite well the propagation of free boundaries. We alsodiscuss the convergence rates of our numerical scheme in the one-component scalar casewhere an explicit solution is known.From a numerical point of view, the multi-component contaminant transport problemhas been studied previously in [1] and [14]. In [1], the authors considered the (non-degenerate) Langmuir isotherm for two-component advection-dominated contaminanttransport using a Local Discontinuous Galerkin Scheme. In [14], the numerical resultswere obtained by adapting the interface-tracking scheme introduced, for scalar porousmedium equation, in [16]. The scheme presented in [16] could be adapted for moregeneral Langmuir and Freundlich type isotherms (including our model problem), if theproblem were first expressed as a system of nonlinear diffusion equations, but the finitespeed of propagation property of free boundaries discussed in [14] can best be capturedif the system has the special diagonal structure as has the problem studied in [14].2. Physical modeling
Let us briefly describe the equations governing the transport of one or more solutesthrough a fluid-saturated porous medium, i.e. a medium characterized by a partitioningof its total volume into a solid phase (solid matrix) and a void or pore space that isfilled by one or more fluids, see [3, 6, 7] for a more detailed description. Assume thatthe flow is at steady state and that the transport is described by advection, moleculardiffusion, mechanical dispersion and chemical reaction (adsorption) between a soluteand the surrounding porous skeleton, and consider the processes on a macroscopic scale.Let Ω ⊂ R d be the domain occupied by the porous medium and let u = u ( x, t ) denotethe solute concentration in the fluid phase. In the one-species model (contamination byone solute) u solves the equation(2.1) ∂ t (cid:0) θu + ρs (cid:1) + ∇ · ( qu − D ∇ u ) = f ( u ) , where θ = θ ( x ) > is the porosity, ρ > is the constant bulk mass density of the solidmatrix, q is the Darcy velocity of the flow ( qu is the advective water flux) and D standsfor the hydrodynamic dispersion matrix including both molecular diffusion and mechan-ical dispersion between the solute and the surrounding porous medium. Moreover, f models source/sink terms that may depend on u and s denotes the concentration of thecontaminant on the porous matrix ( ρ ∂ t s corresponds to the rate of change of concen-tration on the porous matrix due to adsorption or desorption).Normally, s is taken to satisfy an ordinary differential equation of the form(2.2) s t = k F ( u, s ) , where k > is the rate parameter and F is the reaction rate function. In case of fastreaction ( equilibrium adsorption ), we may let k → ∞ in (2.2) and, assuming that the OMPETITIVE CONTAMINANT TRANSPORT 3 equation F ( u, s ) = 0 can be solved for s we obtain s = ψ ( u ) , where ψ = ψ ( u ) is called the adsorption isotherm satisfying, in general, the monotonicityproperty ψ (0) = 0 , ψ is strictly increasing and smooth for u > . Moreover, ψ is said to be of Langmuir type if ψ is strictly concave near u = 0 and lim u → ψ (cid:48) ( u ) < ∞ and of Freundlich type if ψ is strictly concave near u = 0 and lim u → ψ (cid:48) ( u ) = ∞ . The standard Langmuir isotherm is given by ψ ( u ) = N Ku
Ku , K > , where N is the saturation concentration of the adsorbed solute, and the archetypicalFreundlich isotherm is defined by ψ ( u ) = Ku p , K > , with the reaction rate p commonly chosen in (0 , ; the smaller the p , the higher theadsorption at low concentrations is.In the multi-species case (transport of several contaminants), u = ( u , . . . , u m ) is avector-valued function but the evolution of the concentration of each component u i isstill described by (2.1). However, the adsorption process is now competitive (differentspecies competing for the same adsorption sites) thus making the PDE system coupled.The most common multicomponent isotherms are of the form ψ ( u ) = ( ψ ( u ) , ..., ψ m ( u )) , with the components ψ i given by(2.3) ψ i ( u ) = N i K i u i (cid:80) mj =1 K j u j Langmuir , (2.4) ψ i ( u ) = C i (cid:16) m (cid:88) j =1 a ij u j (cid:17) p i − u i Freundlich , see [20, 12, 1, 23, 21, 15]. Above N i , C i , K i , p i and a ij are non-negative constants; N i is the adsorption capacity of the adsorbent i , K i is the Langmuir affinity constant, C i is the Freundlich adsorption constant, p i is the Freundlich adsorption rate and a ij are dimensionless competition coefficients describing the inhibition of species i to theadsorption of species j . By definition a ii = 1 and a ij = 0 if there is no competitionbetween the species i and j .The multicomponent Freundlich isotherm (2.4), first proposed by Sheindorf et al.[20], has been successfully applied to different contaminants, see [22, 12, 23]. However,it is not the only empirical Freundlich type (according to the definition above) isothermthat has been proposed in the literature, see, e.g. , Hinz [13]. Here, we assume that theadsorption process is modelled by such a generalized Freundlich type multicomponentisotherm which can be idealized as ψ i ( u ) = C | u | p − u i , where | · | denotes the usual Euclidean norm in R m . For expediency, we also assumethat the coefficient functions in (2.1) are all constants and that the source term f =( f , . . . , f m ) does not depend on u . Ignoring advection in (2.1), that is setting q = 0 , itfollows that(2.5) θ∂ t (cid:0) u i + ρ C | u | p − u i (cid:1) − D ∆ u i = f i , i = 1 , . . . , m . We may eliminate the constants in (2.5) by redefining the variables through u ∗ i = (cid:16) ρCθ (cid:17) / ( p − u i , t ∗ = tD , x ∗ = θ / xD , f ∗ i = θD (cid:16) ρKθ (cid:17) / ( p − f i . Dropping the stars yields the system of equations(2.6) ∂ t (cid:0) u + | u | p − u (cid:1) − ∆ u = f . System (2.6) is of the form (1.1) , with b ( u ) = u + | u | p − u = ∇ φ ( u ) , where(2.7) φ ( u ) = 12 | u | + 1 p + 1 | u | p +1 , φ : R m → R , BAÍA, BOZORGNIA, MONSAINGEON, AND VIDEMAN is a strictly convex function.3.
Notation and preliminary results
Let us briefly introduce the notation used throughout this work. The Euclidean normin R m is denoted by | · | , C represents a generic positive constant, possibly varying fromline to line, and we often write u ( t )( x ) for u ( t, x ) .Given a real Banach space X , the (Banach) space L p (0 , T ; X ) consists of all measur-able functions u : [0 , T ] → X such that (cid:107) u (cid:107) L p (0 ,T ; X ) = (cid:32)(cid:90) T (cid:107) u ( t ) (cid:107) pX dt (cid:33) p < ∞ , ≤ p < ∞ ,L ∞ (0 , T ; X ) is the space of all measurable u : [0 , T ] → X such that (cid:107) u (cid:107) L ∞ (0 ,T ; X ) = sup t ∈ (0 ,T ) (cid:107) u ( t ) (cid:107) X < ∞ , and the Banach space W ,p (0 , T ; X ) , for ≤ p < ∞ , consists of all u ∈ L p (0 , T ; X ) such that ∂ t u exists in the weak sense and belongs to L p (0 , T ; X ) . We recall that u ∈ L (0 , T ; X ) is weakly differentiable with ∂ t u ∈ L (0 , T ; X ) if and only if for some g ∈ X it holds(3.1) u ( t ) = g + (cid:90) t ∂ t u ( s ) ds a.e. in (0 , T ) . Moreover, if u ∈ W ,p (0 , T ; X ) then u is continuous as a function from [0 , T ] to X , i.e. u ∈ C ([0 , T ]; X ) .We write H k (Ω) = W k, (Ω) and u ∈ H k (Ω; R m ) if each component of u : Ω → R m belongs to H k (Ω) . The action of an element g ∈ H − (Ω; R m ) on ξ ∈ H (Ω; R m ) isdenoted by (cid:104) g, ξ (cid:105) . When no confusion arises, we write (cid:104) g, ξ (cid:105) , instead of (cid:104) g ( t ) , ξ ( t ) (cid:105) , alsofor g ∈ L (0 , T ; H − (Ω; R m )) and ξ ∈ L (0 , T ; H (Ω; R m )) .Given a measurable function u : [0 , T ] → X and h > , we define the backward timedifference quotient of u by ∂ ht u ( t, x ) = u ( t, x ) − u ( t − h, x ) h , ∂ ht u : [ h, T ] → X .
Lemma 3.1.
Let u ∈ W ,p (0 , T ; X ) , ≤ p < ∞ . For any h > , it holds || ∂ ht u || L p ( h,T ; X ) ≤ || ∂ t u || L p (0 ,T ; X ) . Proof.
Given u ∈ W ,p (0 , T ; X ) , it follows from (3.1) that u ( t ) − u ( t − h ) h = 1 h (cid:90) tt − h ∂ t u ( s ) ds, t ∈ [ h, T ] . Hence (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) u ( t ) − u ( t − h ) h (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) X ≤ h (cid:90) tt − h || ∂ t u ( s ) || X ds = 1 h (cid:90) h || ∂ t u ( s + t − h ) || X ds ≤ h q − (cid:16) (cid:90) h || ∂ t u ( s + t − h ) || pX ds (cid:17) p by Hölder’s inequality, with p + q = 1 . Thus, we conclude that (cid:90) Th (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) u ( t ) − u ( t − h ) h (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) pX dt ≤ h (cid:90) h (cid:16) (cid:90) Th || ∂ t u ( s + t − h ) || pX dt (cid:17) ds ≤ h (cid:90) h (cid:16) (cid:90) T || ∂ τ u ( τ ) || pX dτ (cid:17) ds = (cid:90) T || ∂ τ u ( τ ) || pX dτ . using Fubini’s theorem. (cid:3) OMPETITIVE CONTAMINANT TRANSPORT 5
Let h > be a time step. For simplicity take an equidistant partition of the interval [0 T ] and assume that T /h is an integer. Consider a sequence of such time intervals { h k } k ∈ N with h k −→ k →∞ . Again, for simplicity, we ommit the index k and write h → instead of h k −→ k →∞ . Given w ∈ L (0 , T ; X ) , we define w h as the zeroth-order Clémentquasi-interpolant of w . In other words, given n ≡ n ( h ) ∈ { , ..., T /h } we let w ( n ) = 1 h (cid:90) nh ( n − h w ( s ) ds and define w h ∈ L ∞ (0 , T ; X ) through(3.2) w h ( t ) = w ( n ) for t ∈ (( n − h, nh ] , and, if needed, w h ( t ) = w (1) for t ∈ ( − h, . Lemma 3.2.
Given w ∈ L (0 , T ; X ) , then (cid:107) w h (cid:107) L (0 ,T ; X ) ≤ (cid:107) w (cid:107) L (0 ,T ; X ) . Moreover, assuming that ∂ t w ∈ L (0 , T ; X ) ( with obvious modifications if w h is extendedto ( − h, , it holds (3.3) (cid:107) ∂ ht w h (cid:107) L ( h,T ; X ) ≤ (cid:107) ∂ t w (cid:107) L (0 ,T ; X ) . Proof.
From the Jensen’s inequality, it follows that (cid:107) w h (cid:107) L (0 ,T ; X ) = (cid:90) T (cid:107) w h ( t ) (cid:107) X dt = T/h (cid:88) n =1 h (cid:107) w ( n ) (cid:107) X = T/h (cid:88) n =1 h (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) h (cid:90) nh ( n − h w ( s ) ds (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) X ≤ T/h (cid:88) n =1 (cid:90) nh ( n − h || w ( s ) || X ds = (cid:107) w (cid:107) L (0 ,T ; X ) . Similarly, using again the Jensen’s inequality as well as Lemma 3.1, we obtain (cid:107) ∂ ht w h (cid:107) L ( h,T ; X ) = (cid:90) Th (cid:107) ∂ ht w h (cid:107) X dt = T/h (cid:88) n =2 h (cid:13)(cid:13)(cid:13) w ( n ) − w ( n − (cid:13)(cid:13)(cid:13) X = T/h (cid:88) n =2 h (cid:13)(cid:13)(cid:13) h (cid:90) ( n − h ( n − h ( w ( s + h ) − w ( s )) ds (cid:13)(cid:13)(cid:13) X ≤ T/h (cid:88) n =2 (cid:90) ( n − h ( n − h (cid:13)(cid:13)(cid:13)(cid:13) w ( s + h ) − w ( s ) h (cid:13)(cid:13)(cid:13)(cid:13) X ds = (cid:90) T − h (cid:107) ∂ ht w ( s ) (cid:107) X ds = (cid:107) ∂ ht w (cid:107) L (0 ,T − h ; X ) ≤ (cid:107) ∂ t w (cid:107) L (0 ,T ; X ) . (cid:3) We recall (cf. [11]) that a function E : X → [ −∞ , ∞ ] is said to be sequentially lowersemicontinuous, s.l.s.c. for short, if for any { u n } converging (in X ) to u , it holds E ( u ) ≤ lim inf n →∞ E ( u n ) , and that if X is a reflexive Banach space and E ( u ) → ∞ as || u || X → ∞ , then E issequentially coercive with respect to the weak topology of X .Let Ψ : R m → [0 , ∞ ] be the Legendre conjugate of φ , i.e. Ψ( y ) = sup z ∈ R m ( z · y − φ ( z )) and let B : R m → [0 , ∞ ] be given by(3.4) B ( z ) := Ψ( b ( z )) . Given that φ is convex and b ( z ) = ∇ φ ( z ) , we have (see, e.g , [17])(3.5) B ( z ) = z · b ( z ) − φ ( z ) , BAÍA, BOZORGNIA, MONSAINGEON, AND VIDEMAN or, equivalently, since φ (0) = 0 B ( z ) = (cid:90) (cid:16) b ( z ) − b ( sz ) (cid:17) · z ds, z ∈ R m . Lemma 3.3.
Given δ > , it holds (3.6) | b ( z ) | ≤ δB ( z ) + max | w |≤ /δ | b ( w ) | ∀ z ∈ R m . Proof. If | b ( z ) | (cid:54) = 0 , it follows from (3.4) that B ( z ) ≥ b ( z ) · δ | b ( z ) | b ( z ) − φ (cid:18) δ | b ( z ) | b ( z ) (cid:19) ≥ δ | b ( z ) | − max | ζ | =1 /δ φ ( ζ ) ≥ δ | b ( z ) | − δ max | w |≤ /δ | b ( w ) | , where in the last inequality we have used the Mean Value Theorem. If (cid:107) b ( z ) (cid:107) = 0 theresult is trivial. (cid:3) Remark 3.1.
The convexity of φ and (3.5) imply that (3.7) B ( z ) − B ( w ) ≥ (cid:0) b ( z ) − b ( w ) (cid:1) · w ∀ z, w ∈ R m . Problem statement and existence of a weak solution
As mentioned in the Introduction, our aim in this work is to study the solute con-centrations u : [0 , T ] × Ω → R m solving the initial-boundary value problem(4.1) ∂ t b ( u ( t, x )) − ∆ u ( t, x ) = f ( t, x ) in (0 , T ) × Ω u ( t, x ) = u b ( t, x ) on (0 , T ) × Γ b ( u (0 , x )) = b ( x ) in Ω where x = ( x , . . . , x d ) , (cid:52) = (cid:80) di =1 ∂ ∂x i , ∂ t = ∂∂t , Ω ⊂ R d is a smooth domain, Γ = ∂ Ω , < T < ∞ , and f : [0 , T ] × Ω → R m , b : Ω → R m and u b : [0 , T ] × Γ → R m are givenfunctions. We assume that the vector field b : R m → R m can be expressed as a gradientof some non-negative C -convex function φ : R m → R , with b (0) = 0 and φ (0) = 0 , andthat there exists a measurable function u : Ω → R m such that b ( u ) = b . To define a weak solution to problem (4.1) let us assume that u b ∈ L (cid:0) , T ; H / (Γ; R m ) (cid:1) , f ∈ L (cid:0) , T ; H − (Ω; R m ) (cid:1) and b : Ω → R m is such that Ψ( b ) ∈ L (Ω; [0 , ∞ ]) (notethat from Lemma 3.3 this last condition implies that b ∈ L (Ω; R m ) ).Let X be the function space defined through X = (cid:110) v ∈ L (cid:0) , T ; H (Ω; R m ) (cid:1) : b ( v ) ∈ L ∞ (0 , T ; L (Ω; R m )) , ∂ t b ( v ) ∈ L (cid:0) , T ; H − (Ω; R m ) (cid:1)(cid:111) . Definition 4.1.
We say that u ∈ X is a weak solution of (4.1) if a) u | (0 ,T ) × Γ = u b ; b) ∀ ξ ∈ L (0 , T ; H (Ω; R m )) it holds (cid:90) T (cid:10) ∂ t b ( u ) , ξ (cid:11) dt + (cid:90) T (cid:90) Ω ∇ u · ∇ ξ dx dt = (cid:90) T (cid:104) f, ξ (cid:105) dt ; c) ∀ ξ ∈ L (cid:0) , T ; H (Ω; R m ) (cid:1) ∩ W , (cid:0) , T ; L ∞ (Ω; R m ) (cid:1) with ξ ( T ) = 0 , we have (cid:90) T (cid:104) ∂ t b ( u ) , ξ (cid:105) dt + (cid:90) T (cid:90) Ω b ( u ) · ∂ t ξ dx dt + (cid:90) Ω b · ξ (0) dx = 0 . Let us impose the following growth conditions on the vector field b :( H ) There exists C > such that | b ( u ) | ≤ C ( | u | + 1) . ( H ) There exists C > such that | b ( u ) | ≤ C ( B ( u ) + 1) . OMPETITIVE CONTAMINANT TRANSPORT 7
Remark 4.1.
The nonlinear isotherm b ( u ) = u + | u | p − u associated with the competitivecontaminant transport potential (2.7) satisfies the growth conditions ( H ) and ( H ) .Note that in that case B ( u ) = | u | + pp +1 | u | p +1 . We will use Rothe’s method (see [19]) to approximate the initial-boundary value prob-lem (4.1). It means that we will first discretize the time variable and consider a nonlinearelliptic system at each time step. Next, we introduce a convex energy functional andshow that it admits a unique minimizer which also solves the elliptic system. Finally, toconclude that the sequence of semi-discretized solutions converges to a solution of theoriginal problem, we use a compactness argument based on suitable a priori estimatesfor the approximated solutions. The hardest part is to prove compactness in time; herewe follow the approach in [2].Let f h and u h,b be the zeroth-order Clément quasi-interpolants of the functions f and u b , respectively (see (3.2)). Let u (1) ∈ H (Ω; R m ) be such that for all ξ ∈ H (Ω; R m ) itholds(4.2) (cid:90) Ω b ( u (1) ( x )) − b ( x ) h · ξ ( x ) dx + (cid:90) Ω ∇ u (1) ( x ) · ∇ ξ ( x ) dx = (cid:104) f (1) , ξ (cid:105) and u (1) | Γ = u (1) b . In general, for n = 2 , ..., T /h , and given u ( n − ∈ H (Ω; R m ) , let u ( n ) ∈ H (Ω; R m ) be a solution to the problem(4.3) (cid:90) Ω b ( u ( n ) ( x )) − b ( u ( n − ( x )) h · ξ ( x ) dx + (cid:90) Ω ∇ u ( n ) ( x ) · ∇ ξ ( x ) dx = (cid:104) f ( n ) , ξ (cid:105) for all ξ ∈ H (Ω; R m ) , with u ( n ) | Γ = u ( n ) b . The existence of such functions u ( n ) isestablished in the following lemma (see also Remark 4.2). Lemma 4.1.
Let h > and assume that u ∈ L (Ω; R m ) , g ∈ H − (Ω; R m ) , w ∈ H / (Γ; R m ) and that hypothesis ( H ) holds. Then the variational problem: find u ∈ H (Ω; R m ) such that u | Γ = w and (4.4) (cid:90) Ω b ( u ( x )) − b ( u ( x )) h · ξ ( x ) dx + (cid:90) Ω ∇ u ( x ) · ∇ ξ ( x ) dx = (cid:104) g, ξ (cid:105) , for all ξ ∈ H (Ω; R m ) , admits at least one solution.Proof. Let E : H (Ω; R m ) → [ −∞ , ∞ ] be a functional defined through E ( v ) = (cid:90) Ω φ ( v ( x ) + w ( x )) − b ( u ( x )) · ( v ( x ) + w ( x )) h dx + 12 (cid:90) Ω |∇ ( v ( x ) + w ( x )) | dx − (cid:104) g, v + w (cid:105) , where w ∈ H (Ω; R m ) is such that w coincides with w in the trace sense on Γ . Observethat hypothesis ( H ) , the fact that b = ∇ φ , and the assumption on u guarantee that | E ( v ) | < ∞ for all v ∈ H (Ω; R m ) . Moreover, if v is a smooth critical point of E then u = v + w satisfies the Euler-Lagrange equation (4.4) with u | Γ = w .Recalling that φ is a non-negative function, the Poincaré, Cauchy-Schwarz and Younginequalities yield E ( v ) ≥ C (cid:104) (cid:107) v (cid:107) H − (cid:107) v (cid:107) H − (cid:105) , for some C = C ( h, (cid:107) u (cid:107) L , (cid:107) w (cid:107) H , (cid:107) g (cid:107) H − ) which shows that E is coercive and thusminimizing sequences { v k } are bounded in H . In addition, the functional E is convexand lower semicontinuous, thus weakly lower semicontinuous, and H (Ω) is a reflexiveBanach space so there exists a subsequence converging to an element v (weakly in H and strongly in L ) which minimizes E . Moreover, u = ( v + w ) ∈ H (Ω; R m ) is asolution to (4.4) with u | Γ = w . (cid:3) Remark 4.2.
The existence of a solution u ( n ) to problem (4.3) , for n ∈ { , ..., T /h } ,follows from Lemma 4.1 with g = f ( n ) and w = u ( n ) b , while the existence of a solution BAÍA, BOZORGNIA, MONSAINGEON, AND VIDEMAN u (1) to problem (4.2) can be established by replacing hypothesis ( H ) by ( H ) in theprevious lemma. Note that ( H ) implies b ∈ L (Ω; R m ) . Remark 4.3.
The solution in Lemma 4.1 is unique since φ is strictly convex and thus E is strictly convex. Let us now construct the approximated solution of (4.1) globally for all t ∈ [0 , T ] .Let u h ∈ L ∞ (0 , T ; H (Ω)) be the (piecewise constant in time) function defined by(4.5) u h ( t, x ) := u ( n ) ( x ) for ( t, x ) ∈ (( n − h, nh ] × Ω , n ∈ { , ..., T /h } , where u ( n ) is the solution to problem (4.2) if n = 1 or to problem (4.3) for n ∈{ , ..., T /h } . We extend u h ( t ) to ( − h, by u .By construction, u h satisfies the following result. Lemma 4.2.
Let h > and assume that hypotheses ( H ) and ( H ) are valid. Then,for all t ∈ (0 , T ) and all ξ ∈ H (Ω) , the function u h satisfies (4.6) (cid:90) Ω ∂ ht b ( u h ( t )) · ξ dx + (cid:90) Ω ∇ u h ( t ) · ∇ ξ dx = (cid:104) f h ( t ) , ξ (cid:105) . Moreover, u h | (0 ,T ) × Γ = u h,b and for all ξ ∈ L (cid:0) , T ; H (Ω; R m ) (cid:1) it holds (cid:90) T (cid:90) Ω ∂ ht b ( u h ( t )) · ξ ( t ) dx dt + (cid:90) T (cid:90) Ω ∇ u h ( t ) · ∇ ξ ( t ) dx dt = (cid:90) T (cid:104) f h ( t ) , ξ ( t ) (cid:105) dt . We also have the following "integration by parts" formula.
Lemma 4.3.
Let u h ∈ L ∞ (cid:0) , T ; H (Ω) (cid:1) and assume that conditions ( H ) and ( H ) arevalid. Then (cid:90) T (cid:90) Ω ∂ ht b ( u h ( t )) · ξ ( t ) dx dt = − (cid:90) T − h (cid:90) Ω b ( u h ( t )) · ∂ ht ξ ( t + h ) dx dt + 1 h (cid:32)(cid:90) TT − h (cid:90) Ω b ( u h ( t )) · ξ ( t ) dx dt − (cid:90) h (cid:90) Ω b · ξ ( t ) dx dt (cid:33) , for all ξ ∈ L (cid:16) , T ; H (Ω; R m ) (cid:17) ∩ W , (cid:16) , T ; L ∞ (Ω; R m ) (cid:17) .Proof. It is easy to see that (cid:90) T (cid:90) Ω ∂ ht b ( u h ( t )) · ξ ( t ) dx dt = T/h (cid:88) n =1 (cid:90) nh ( n − h (cid:90) Ω b ( u ( n ) ) − b ( u ( n − ) h · ξ ( t ) dx dt = − T/h − (cid:88) n =1 (cid:90) nh ( n − h (cid:90) Ω b ( u ( n ) ) · ξ ( t + h ) − ξ ( t ) h dx dt + 1 h (cid:18) (cid:90) TT − h (cid:90) Ω b ( u ( n ) ) · ξ ( t ) dx dt − (cid:90) h (cid:90) Ω b ( u (0) ) · ξ ( t ) dx dt (cid:19) = − (cid:90) T − h (cid:90) Ω b ( u h ( t ) · ∂ ht ξ ( t + h ) dx dt + 1 h (cid:18) (cid:90) TT − h (cid:90) Ω b ( u h ( t )) · ξ ( t ) dx dt − (cid:90) h (cid:90) Ω b · ξ ( t ) dx dt (cid:19) . (cid:3) We will now present a series of auxiliary lemmas needed in proving the main existenceresult. We will start by establishing an a priori estimate for smooth solutions. For thatwe need the following regularity assumption on the boundary data.
Assumption 4.1.
Assume that the extension of u b ∈ L (0 , T ; H / (Ω; R m )) to Ω , stilldenoted by u b , is such that u b ∈ L ∞ (0 , T ; H (Ω; R m )) and ∂ t u b ∈ L (0 , T ; L (Ω; R m )) . OMPETITIVE CONTAMINANT TRANSPORT 9
Lemma 4.4.
Assume that problem (4.1) admits a smooth solution and that Assumption4.1 is valid. Then the following a priori estimate holds (4.7) (cid:107) B ( u ) (cid:107) L ∞ (0 ,T ; L ) + (cid:107) u (cid:107) L (0 ,T ; H ) ≤ C (cid:0) (cid:107) Ψ( b ) (cid:107) L + (cid:107) f (cid:107) L (0 ,T ; H − ) + (cid:107) u b (cid:107) L (0 ,T ; H ) + (cid:107) ∂ t u b (cid:107) L (0 ,T ; L ) ) , for some C = C (Ω , T ) .Proof. Multiplying the first equation in (4.1) by u − u b and integrating over Ω , gives(4.8) (cid:90) Ω ∂ t b ( u ) · ( u − u b ) dx + (cid:90) Ω |∇ u | dx = (cid:90) Ω ∇ u · ∇ u b dx + (cid:104) f, u − u b (cid:105) . On the other hand, from (3.5) it follows that(4.9) ∂ t B ( u ) = ∂ t b ( u ) · u + b · ∂ t u − ∇ φ · ∂ t u = ∂ t b ( u ) · u . Substituting (4.9) into (4.8), yields(4.10) (cid:90) Ω ∂ t B ( u ) dx + (cid:90) Ω |∇ u | dx = (cid:90) Ω ∂ t b ( u ) · u b dx + (cid:90) Ω ∇ u · ∇ u b dx + (cid:104) f, u − u b (cid:105) . Integrating (4.10) over (0 , t ) , recalling (3.4) and using the initial condition b (0 , x ) = b ( x ) , we obtain(4.11) (cid:90) Ω B ( u ( t )) dx + (cid:90) t (cid:90) Ω |∇ u | dx ds = (cid:90) Ω Ψ( b ) dx + (cid:90) t (cid:90) Ω ∂ t b ( u ) · u b dx ds + (cid:90) t (cid:90) Ω ∇ u · ∇ u b dx ds + (cid:90) t (cid:104) f, u − u b (cid:105) ds . Using Hölder’s, Young’s and Poincaré’s inequalities, one easily shows that(4.12) (cid:90) t (cid:90) Ω ∇ u · ∇ u b dx ds + (cid:90) t (cid:104) f, u − u b (cid:105) ds ≤ C (cid:16) (cid:107) f (cid:107) L (0 ,T ; H − ) + (cid:107) u b (cid:107) L (0 ,T ; H ) (cid:17) + (cid:15) (cid:90) t (cid:90) Ω |∇ u | dx ds . holds for all t ∈ [0 , T ] and for any (cid:15) > and where C = C (Ω , (cid:15) − ) . Moreover, integrationby parts gives (cid:90) t (cid:90) Ω ∂ t b ( u ) · u b dx ds = (cid:90) Ω [ b ( u ( t )) u b ( t ) − b ( u (0)) u b (0)] dx − (cid:90) t (cid:90) Ω b ( u ) · ∂ t u b dx ds, from where, using again Hölder’s and Young’s inequalities as well as hypothesis ( H ) , itfollows, for all t ∈ [0 , T ] and for any (cid:15) > , that(4.13) (cid:90) t (cid:90) Ω ∂ t b ( u ) · u b dx ds ≤ C (cid:16) (cid:107) Ψ( b ) (cid:107) L + (cid:107) u b (cid:107) L ∞ (0 ,T : L ) + (cid:107) ∂ t u b (cid:107) L (0 ,T ; L ) + (cid:15) max t ∈ [0 ,T ] (cid:90) Ω B ( u ( t )) dx (cid:17) , where C = C (Ω , (cid:15) − , T ) . Using estimates (4.12) and (4.13) in (4.11) and choosing (cid:15) > small enough, we first obtain(4.14) (cid:90) Ω B ( u ( t )) dx + (cid:90) t (cid:90) Ω |∇ u | dx ds ≤ C (cid:16) (cid:107) Ψ( b ) (cid:107) L + (cid:107) f (cid:107) L (0 ,T ; H − ) + (cid:107) u b (cid:107) L (0 ,T ; H ) + (cid:107) ∂ t u b (cid:107) L (0 ,T ; L ) + (cid:15) max t ∈ [0 ,T ] (cid:90) Ω B ( u ( t )) dx (cid:17) . Finally, taking the supremum over t ∈ [0 , T ] in (4.14), choosing (cid:15) > small enough andusing the Poincaré’s inequality for u − u b , yields estimate (4.7). (cid:3) The next lemmas concern the semi-discrete solution u h defined in (4.5). We will firstprove an a priori bound, uniform in h > , similar to the one shown in Lemma 4.4. Lemma 4.5.
Assume that hypotheses ( H ) and ( H ) are valid and Assumption 4.1holds. There exists C = C (Ω) such that for any h > (4.15) sup t ∈ [0 ,T ] (cid:90) Ω B ( u h ( t )) dx + (cid:90) T (cid:90) Ω |∇ u h | dx dt ≤ C (cid:0) (cid:107) Ψ( b ) (cid:107) L + (cid:107) f (cid:107) L (0 ,T ; H − ) + (cid:107) u b (cid:107) L (0 ,T ; H ) + (cid:107) ∂ t u b (cid:107) L (0 ,T ; L ) ) . Proof.
Given any n ∈ { , ..., T /h } , consider the difference ξ = u ( n ) − u ( n ) b ∈ H (Ω) as atest function in (4.3) (or in (4.2)). In view of (3.7), we obtain for all n ≥ (cid:90) Ω (cid:0) B ( u ( n ) ) − B ( u ( n − ) (cid:1) dx + h (cid:90) Ω |∇ u ( n ) | dx ≤ h (cid:90) Ω ∇ u ( n ) · ∇ u ( n ) b dx + h (cid:68) f ( n ) , u ( n ) − u ( n ) b (cid:69) + (cid:90) Ω (cid:16) b ( u ( n ) ) − b ( u ( n − ) (cid:17) · u ( n ) b dx. Using Poincaré’s and Young’s inequalities, we immediately see that(4.16) (cid:90) Ω (cid:0) B ( u ( n ) ) − B ( u ( n − ) (cid:1) dx + h (cid:90) Ω |∇ u ( n ) | dx ≤ Ch (cid:16) (cid:107)∇ u ( n ) b (cid:107) L + (cid:107) f ( n ) (cid:107) H − (cid:17) + (cid:90) Ω (cid:16) b ( u ( n ) ) − b ( u ( n − ) (cid:17) · u ( n ) b dx . Adding inequalities (4.16) from n = 1 to n = N , letting t = N h , with N ∈ { , ..., T /h } ,and recalling that u h is piecewise constant in time, gives(4.17) (cid:90) Ω B ( u h ( t )) dx + 12 (cid:90) t (cid:90) Ω |∇ u h ( t ) | dx ≤ (cid:90) Ω B ( u ) dx + C (cid:16) (cid:107)∇ u b (cid:107) L (0 ,T ; L ) + (cid:107) f (cid:107) L (0 ,T ; H − ) (cid:17) + N (cid:88) n =1 (cid:90) Ω (cid:16) b ( u ( n ) ) − b ( u ( n − ) (cid:17) · u ( n ) b dx . where we have also taken Lemma 3.2 into account.Since for any pair of sequences { x n } n ∈ N and { y n } n ∈ N the “summation by parts”formula N (cid:88) n =1 ( x n − x n − ) y n = ( x N y N − x y ) − h N (cid:88) n =1 x n − (cid:18) y n − y n − h (cid:19) holds, we can write S := N (cid:88) n =1 (cid:90) Ω (cid:16) b ( u ( n ) ) − b ( u ( n − ) (cid:17) · u ( n ) b dx = (cid:90) Ω b ( u h ( t )) · u h,b ( t ) dx − (cid:90) Ω b ( u ) · u h,b (0) dx − t (cid:90) (cid:90) Ω b ( u h ( s − h, x )) ∂ ht u h,b ( s, x ) dx ds. Hence, using Hölder’s and Young’s inequalities, Lemma 3.2 and the fact that u h ispiecewise constant in time, we get(4.18) | S | ≤ ε (cid:16) (cid:107) b ( u h ( t )) (cid:107) L + (cid:107) b ( u h ) (cid:107) L ( − h,t : L ) (cid:17) + C (cid:16) (cid:107) u b (cid:107) L ∞ (0 ,T ; L ) + (cid:107) ∂ ht u b (cid:107) L (0 ,T ; L ) (cid:17) + (cid:90) Ω | b ( u ) | dx . Moreover, in view of hypotheses ( H ) and ( H ) (4.19) (cid:107) b ( u h ( t )) (cid:107) L ≤ C (cid:18) (cid:90) Ω B ( u h ( t )) dx (cid:19) OMPETITIVE CONTAMINANT TRANSPORT 11 and(4.20) (cid:107) b ( u h ) (cid:107) L (( − h,t ) × Ω) = h (cid:107) b ( u ) (cid:107) L + (cid:107) b ( u h ) (cid:107) L (0 ,t,L ) ≤ C (cid:16) (cid:107) B ( u ) (cid:107) L + (cid:107) u h (cid:107) L (0 ,t ; L ) (cid:17) ≤ C (cid:16) (cid:107) Ψ( u ) (cid:107) L + (cid:107) u b (cid:107) L (0 ,T ; H ) + (cid:107)∇ u h (cid:107) L (0 ,t ; L ) (cid:17) , where we have also use Poincaré’s inequality for u h − u h,b ∈ H (Ω) and Lemma 3.2 .Collecting estimates (4.18), (4.19), (4.20) and using them in (4.17) and choosing (cid:15) > small enough, gives finally (cid:90) Ω B ( u h ( t )) dx + (cid:90) t (cid:90) Ω |∇ u h | dx dt ≤ C (cid:0) (cid:107) Ψ( b ) (cid:107) L + (cid:107) f (cid:107) L (0 ,T ; H − ) + (cid:107) u b (cid:107) L (0 ,T ; H ) + (cid:107) ∂ t u b (cid:107) L (0 ,T ; L ) ) , with some constant C = C (Ω) > , independent of h , from which estimate (4.15) readilyfollows. (cid:3) Remark 4.4.
Using Poincaré’s inequality for and Lemma 3.2, we also obtain (cid:107) u h (cid:107) L (0 ,T ; H ) ≤ C (cid:0) (cid:107) Ψ( b ) (cid:107) L + (cid:107) f (cid:107) L (0 ,T ; H − ) + (cid:107) u b (cid:107) L (0 ,T ; H ) + (cid:107) ∂ t u b (cid:107) L (0 ,T ; L ) ) . The following lemma is a straightforward modification of Lemma 1.9 presented in [2].It will be essential in proving the main existence theorem. Its proof, which we omit, isbased on Minty’s monotonicity argument and makes use of estimate (3.6) of Lemma 3.3.
Lemma 4.6.
Assume that { u h } converges weakly in L (0 , T ; H (Ω)) to some u ∈ L (0 , T ; H (Ω)) and let τ ∈ (0 , T ) . Moreover, assume that there exists a constant C > such that T − τ (cid:90) (cid:90) Ω (cid:0) b ( u h ( s + τ )) − b ( u h ( s )) (cid:1) · (cid:0) u h ( s + τ ) − u h ( s ) (cid:1) dx ds ≤ C max { τ, √ τ } and that (cid:90) Ω B ( u h ( t )) dx ≤ C , for all < t < T . Then b ( u h ) converges to b ( u ) in L (0 , T ; L (Ω)) and B ( u h ) converges to B ( u ) almosteverywhere. Lemma 4.5 guarantees that u h satisfies the second assumption of Lemma 4.6. Thefirst assumption is established in the next two lemmas. Lemma 4.7.
Let h > and t , t ∈ (0 , T ) be given. Then for all ξ ∈ H (Ω) it holds (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:90) Ω [ b ( u h ( t )) − b ( u h ( t )] · ξ dx + t (cid:90) t (cid:18)(cid:90) Ω ∇ u h ( s ) · ∇ ξ dx − (cid:104) f h ( s ) , ξ (cid:105) (cid:19) ds (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ √ h (cid:0) (cid:107)∇ u h (cid:107) L (0 ,T ; L ) + (cid:107) f (cid:107) L (0 ,T ; H − ) (cid:1) (cid:107)∇ ξ (cid:107) L . Proof.
From (4.3)-(4.2) it follows that(4.21) (cid:90) Ω (cid:104) b (cid:16) u ( n ) (cid:17) − b (cid:16) u ( n − (cid:17)(cid:105) · ξ dx + h (cid:90) Ω ∇ u ( n ) · ∇ ξ dx = h (cid:68) f ( n ) , ξ (cid:69) . for all n ≥ . Let N , N ∈ { , ..., T /h } be such that t ∈ (( N − h, N h ] and t ∈ (( N − h, N h ] . Without loss of generality assume that N < N . Summing from N + 1 to N in (4.21), gives N (cid:88) n = N +1 (cid:90) Ω (cid:104) b (cid:16) u ( n ) (cid:17) − b (cid:16) u ( n − (cid:17)(cid:105) · ξ dx + h N (cid:88) n = N +1 (cid:90) Ω ∇ u ( n ) · ∇ ξ dx = h N (cid:88) n = N +1 (cid:104) f n , ξ (cid:105) , In view of the definition of u h and f h and observing that u h ( t i ) = u ( N i ) , i = 1 , , weobtain(4.22) (cid:90) Ω [ b ( u h ( t )) − b ( u h ( t )) · ξ dx = t (cid:90) t (cid:90) Ω −∇ u h ( s ) · ∇ ξ dx + (cid:104) f h ( s ) , ξ (cid:105) ds − N h (cid:90) t (cid:90) Ω −∇ u h ( s ) · ∇ ξ dx + (cid:104) f h ( s ) , ξ (cid:105) ds + N h (cid:90) t (cid:90) Ω −∇ u h ( s ) · ∇ ξ dx + (cid:104) f h ( s ) , ξ (cid:105) ds. Using Hölder’s inequality and Lemma 3.2, we can estimate(4.23) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N i h (cid:90) t i (cid:90) Ω −∇ u h ( s ) · ∇ ξ dx + (cid:104) f h ( s ) , ξ (cid:105) ds (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = N i h (cid:90) t i (cid:107)∇ u h ( s ) (cid:107) L (cid:107)∇ ξ (cid:107) L ds + N i h (cid:90) t i (cid:107) f h ( s ) (cid:107) H − (cid:107)∇ ξ (cid:107) L ds ≤ √ h (cid:0) (cid:107)∇ u h (cid:107) L (0 ,T ; L ) + (cid:107) f (cid:107) L (0 ,T ; H − ) (cid:1) (cid:107)∇ ξ (cid:107) L , i = 1 , , and the result follows from (4.22) and (4.23). (cid:3) Lemma 4.8.
Given τ ∈ (0 , T ) there exists a constant C = C (Ω , T, f, u , u b , τ ) such thatfor all h > it holds T − τ (cid:90) (cid:90) Ω [ b ( u h ( t + τ )) − b ( u h ( t ))] · [ u h ( t + τ ) − u h ( t )] dx dt ≤ C max { τ, √ τ } . Proof.
Consider first τ < h . Since u h ( t + τ ) ≡ u h ( t ) ∀ t ∈ (( n − h, nh − τ ] , we get T − τ (cid:90) (cid:90) Ω [ b ( u h ( t + τ )) − b ( u h ( t ))] · [ u h ( t + τ ) − u h ( t )] dx dt = T/h − (cid:88) n =1 nh (cid:90) nh − τ (cid:90) Ω [ b ( u h ( t + τ )) − b ( u h ( t ))] · [ u h ( t + τ ) − u h ( t )] dx dt = T/h − (cid:88) n =1 nh (cid:90) nh − τ (cid:90) Ω (cid:104) b ( u ( n +1) ) − b ( u ( n ) ) (cid:105) · (cid:104) u ( n +1) − u ( n ) (cid:105) dx dt = τ T/h − (cid:88) n =1 (cid:90) Ω (cid:104) b ( u ( n +1) ) − b ( u ( n ) ) (cid:105) · (cid:104) u ( n +1) − u ( n ) (cid:105) dx = τ T/h − (cid:88) n =1 (cid:90) Ω (cid:104) b ( u ( n +1) ) − b ( u ( n ) ) (cid:105) · (cid:104) v ( n +1) − v ( n ) (cid:105) dx − τ T/h − (cid:88) n =1 (cid:90) Ω (cid:104) b ( u ( n +1) ) − b ( u ( n ) ) (cid:105) · (cid:104) u ( n +1) b − u ( n ) b (cid:105) dx := τ ( I ,h + I ,h ) , OMPETITIVE CONTAMINANT TRANSPORT 13 where v ( n ) = u ( n ) − u ( n ) b ∈ H (Ω) . Equation (4.6) implies that I ,h = T/h − (cid:88) n =1 h (cid:104) (cid:104) f ( n ) , v ( n +1) − v ( n ) (cid:105) − (cid:90) Ω ∇ u n +1 · ∇ ( v ( n +1) − v ( n ) ) dx (cid:105) . From Lemma 3.2 and estimate (4.15) it then follows that I .h ≤ C (cid:0) (cid:90) T (cid:107) f h (cid:107) H − (cid:107)∇ v h (cid:107) L dt + (cid:90) T (cid:107)∇ u h (cid:107) L (cid:107)∇ v h (cid:107) L dt (cid:1) ≤ K , for some constant K = K (Ω , f, b , u b ) > , independent of h and τ . Moreover, I ,h = 1 h (cid:90) T − h (cid:90) Ω [ b ( u h ( t + h )) − b ( u h ( t ))] · [ u h,b ( t + h ) − u h,b ( t )] dx dt ≤ C (cid:107) b ( u h ) (cid:107) L (0 ,T ; L ) (cid:107) ∂ ht u h,b (cid:107) L (0 ,T ; L ) ≤ C (1 + (cid:107) u h (cid:107) L (0 ,T ; L ) ) (cid:107) ∂ t u b (cid:107) L (0 ,T ; L ) ≤ K , with some K = K (Ω , T, f, b , u b ) and where we have taken account hypothesis ( H ) as well as estimates (3.3) and (4.15). We thus conclude that T − τ (cid:90) (cid:90) Ω [ b ( u h ( t + τ )) − b ( u h ( t ))] · [ u h ( t + τ ) − u h ( t )] dx dt ≤ C τ, for some K = K (Ω , T, f, b , u b , τ ) .Assume next that τ ≥ h and fix t ∈ [0 , T − τ ] . In view of Lemma 4.7 and estimate(4.15), for all ξ ( t ) ∈ H (Ω) it holds (cid:90) Ω (cid:104) b ( u h ( t + τ ) − b ( u h ( t )) (cid:105) · ξ ( t ) dx + t + τ (cid:90) t (cid:16) (cid:90) Ω ∇ u h ( s ) · ∇ ξ dx − (cid:104) f h ( s ) , ξ ( t ) (cid:105) ds (cid:17) ≤ C √ h (cid:107)∇ ξ ( t ) (cid:107) L ≤ C √ τ (cid:107)∇ ξ ( t ) (cid:107) L , where the constant C > is independent of τ and h . Taking ξ ( t ) = v h ( t + τ ) − v h ( t ) ,where v h ( t ) = u h ( t ) − u h,b ( t ) , and integrating in t from to T − τ , we get T − τ (cid:90) (cid:90) Ω ( b ( u h ( t + τ )) − b ( u h ( t ))) · ( u h ( t + τ ) − u h ( t )) dx dt ≤ C (cid:112) τ ( T − τ ) + T − τ (cid:90) (cid:90) Ω ( b ( u h ( t + τ )) − b ( u h ( t ))) · ( u h,b ( t + τ ) − u h,b ( t )) dx dt − T − τ (cid:90) (cid:16) t + τ (cid:90) t (cid:90) Ω ∇ u h ( s ) · ∇ ξ ( t ) dx ds (cid:17) dt + T − τ (cid:90) (cid:16) t + τ (cid:90) t (cid:104) f h ( s ) , ξ ( t ) (cid:105) ds (cid:17) dt := C (cid:112) τ ( T − τ ) + I ,h + I ,h + I ,h , where we have also used estimate (4.15) to bound (cid:107)∇ ξ (cid:107) L (0 ,T ; L ) .The integral I ,h can be estimated as I ,h in the previous step yielding | I ,h | ≤ C τ .
As for I ,h , we proceed as follows |I ,h | ≤ (cid:90) T − τ (cid:16) (cid:107)∇ ξ ( t ) (cid:107) L (cid:90) t + τt (cid:107)∇ u h ( s ) (cid:107) L ds (cid:17) dt ≤ √ τ (cid:107)∇ u h (cid:107) L (0 ,T ; L ) (cid:90) T − τ (cid:107)∇ ξ ( t ) (cid:107) L dt ≤ C (cid:112) τ ( T − τ ) , where we have used again estimate (4.15). For the integral I ,h we obtain similarly |I ,h | ≤ C (cid:112) τ ( T − τ ) , and we thus conclude that T − τ (cid:90) (cid:90) Ω ( b ( u h ( t + τ )) − b ( u h ( t ))) · ( u h ( t + τ ) − u h ( t )) dx dt ≤ C (cid:0) τ + (cid:112) τ ( T − τ ) (cid:1) ≤ C max { τ, √ τ } , for some constant C = C (Ω , T, f, b , u b ) independent of h and τ . (cid:3) We are now in a position to prove the main existence theorem.
Theorem 1.
Assume that hypotheses ( H ) and ( H ) are valid and that Assumption 4.1holds. There exists a unique weak solution u ∈ X to problem (4.1) such that B ( u ) ∈ L ∞ (0 , T ; L (Ω; R m )) .Proof. Let us start with existence: for any h > , let u h ∈ L (0 , T ; H (Ω)) be thepiecewise constant function defined in (4.5) and which, in view of Lemmas 4.2 and 4.3,satisfiesi) u h | (0 ,T ) × Γ = u h,b ;ii) for all ξ ∈ L (cid:0) , T ; H (Ω; R m ) (cid:1) (4.24) (cid:90) T (cid:90) Ω ∂ ht b ( u h ( t )) · ξ ( t ) dx dt + (cid:90) T (cid:90) Ω ∇ u h ( t ) · ∇ ξ ( t ) dx dt = (cid:90) T (cid:104) f h ( t ) , ξ ( t ) (cid:105) dt ; iii) for all ξ ∈ L (cid:0) , T ; H (Ω; R m ) (cid:1) ∩ W , (cid:0) , T ; L ∞ (Ω; R m ) (cid:1) , with ξ ( T ) = 0 ,(4.25) (cid:90) T (cid:90) Ω ∂ ht b ( u h ( t )) · ξ ( t ) dx dt = − (cid:90) T − h (cid:90) Ω b ( u h ( t )) ∂ ht ξ ( t + h ) dx dt + 1 h (cid:34)(cid:90) TT − h (cid:90) Ω b ( u h ( t )) · ξ ( t ) dx dt − (cid:90) − h (cid:90) Ω b · ξ ( t + h ) dx dt (cid:35) . Consider a sequence { u h } of such functions and note that Lemma 4.5 (see also Remark4.4) guarantees the existence of a function u ∈ L (0 , T ; H (Ω; R m )) such that (up tosubsequence)(4.26) u h (cid:42) h → u in L (0 , T ; H (Ω; R m )) . For u to be a weak solution to problem (4.1), it needs to meet the conditions a), b) and c) of Definition 4.1 and be such that b ( u ) ∈ L ∞ (0 , T ; L (Ω; R m )) and ∂ t b ( u ) ∈ L (0 , T ; H − (Ω; R m ) . Now, u h,b is the (piecewise constant in time) Clément interpolant of u b and u h satisfies(4.26) so u h,b − u h (cid:42) h → u b − u in L (0 , T ; H (Ω , R m )) . Hence, u b − u ∈ L (0 , T ; H (Ω , R m )) and condition a holds.To show that b holds, we first observe that f h being the Clément interpolant of f and (4.26) guarantee that (cid:90) T (cid:104) f h , ξ (cid:105) dt −→ h → (cid:90) T (cid:104) f, ξ (cid:105) dt , (cid:90) T ∇ u h · ξ dt −→ h → (cid:90) T ∇ u · ξ dt for all ξ ∈ L (0 , T ; H (Ω; R m )) . Moreover, from (4.24), we obtain (cid:12)(cid:12)(cid:12) (cid:90) T (cid:10) ∂ ht b ( u h ) , ξ (cid:11) dt (cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12) (cid:90) T (cid:90) Ω ∂ ht b ( u h ( t )) · ξ ( t ) dx dt (cid:12)(cid:12)(cid:12) ≤ ( (cid:107)∇ u h (cid:107) L (0 ,T ; L ) + (cid:107) f h (cid:107) L (0 ,T,H − ) ) (cid:107) ξ (cid:107) L (0 ,T ; H ) , for all ξ ∈ L (0 , T ; H (Ω; R m )) . Using Lemmas 4.5 and 3.2, we thus conclude that(4.27) (cid:107) ∂ ht b ( u h ) || L (0 ,T ; H − ) ≤ C , for some constant
C > depending on f, u b , b and Ω but independent of h . OMPETITIVE CONTAMINANT TRANSPORT 15
Estimate (4.27) guarantees the existence of λ ∈ L (0 , T ; H − (Ω; R m )) such that ∂ ht b ( u h ) (cid:42) λ in L (0 , T ; H − (Ω; R m )) . To conclude b) , it remains to be checked that λ = ∂ t b ( u ) in L (0 , T ; H − (Ω; R m )) . Tothis end, it is enough to show that we can pass to the limit in (4.25). The left-hand sideimmediately yields (cid:90) T (cid:90) Ω ∂ ht b ( u h ( t )) · ξ ( t ) dx dt = (cid:90) T (cid:10) ∂ ht b ( u h ) , ξ (cid:11) −→ h → (cid:90) T (cid:104) λ, ξ (cid:105) For the first term on the right-hand side, the estimate (cid:107) u h (cid:107) L (0 ,T ; H ) ≤ C from Re-mark 4.4 allows to apply Lemma 4.6 and conclude that b ( u h ) → b ( u ) strongly in L (0 , T ; L (Ω; R m )) and, in particular, b ( u h ) → b ( u ) pointwise a.e. in (0 , T ) × Ω . Since,obviously ∂ ht ξ ( t + h, x ) = ξ ( t + h,x ) − ξ ( t,x ) h → ∂ t ξ ( t, x ) a.e. in (0 , T ) × Ω , and by Lemma 3.1 (cid:107) ∂ ht ξ (cid:107) L (0 ,T ; L ∞ ) ≤ C , we have (cid:90) T − h (cid:90) Ω b ( u h ( t )) ∂ ht ξ ( t + h ) dx dt → (cid:90) T (cid:90) Ω b ( u ( t )) ∂ t ξ ( t + h ) dx dt . For the second term on the right-hand side in (4.25), the regularity ξ ∈ W , (0 , T ; L ∞ (Ω; R m ) guarantees that h (cid:90) TT − h (cid:107) ξ ( t ) (cid:107) L ∞ (Ω) dt → (cid:107) ξ ( T ) (cid:107) L ∞ (Ω) = 0 , thus (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) h (cid:90) TT − h (cid:90) Ω b ( u h ( t )) · ξ ( t ) dx dt (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:107) b ( u h ) (cid:107) L ∞ (0 ,T ; L ) h (cid:90) TT − h (cid:107) ξ ( t ) (cid:107) L ∞ (Ω) dt → . Moreover, h (cid:82) h ξ ( t ) dt → ξ (0) strongly in L ∞ (Ω) , and thus h (cid:90) h (cid:90) Ω b · ξ ( t ) dx dt = (cid:90) Ω b (cid:32) h (cid:90) h ξ ( t ) dt (cid:33) dx → (cid:90) Ω b ξ (0) dx. Therefore, passing to the limit in (4.25), we conclude that (cid:90) T (cid:104) λ, ξ ( t ) (cid:105) dt = − (cid:90) T (cid:90) Ω b ( u ( t )) ∂ t ξ ( t ) dx dt − (cid:90) Ω b · ξ (0) dx, for all ξ ∈ L (cid:16) , T ; H (Ω; R m ) (cid:17) ∩ W , (cid:16) , T ; L ∞ (Ω; R m ) (cid:17) with ξ ( T ) = 0 , which meansof course that λ = ∂ t b ( u ) ∈ L (0 , T ; H − (Ω; R m ) and the initial condition b ( u (0)) = b are satisfied in the weak sense. Hence, conditions b and c hold.As for the uniqueness, the proof of [2, Theorem 2.4] applies directly since we areconsidering linear diffusion in (1.1); note that hypothesis ( H ) guarantees that b ( u ) ∈ L (0 , T ; L (Ω; R m )) .Finally, we show that B ( u ) ∈ L ∞ (0 , T ; L (0 , T ; R m )) from which, using hypotheses ( H ) , it follows that b ( u ) ∈ L ∞ (0 , T ; L (Ω; R m )) so that u ∈ X . (cid:3) Positive properties of the solutions
Under an additional structural assumption on the contaminant transport potential φ , we retrieve nonnegative solutions provided the source data is (componentwise) non-negative. More precisely, we can state the following result. Theorem 2.
Let the assumptions of Theorem 1 hold and assume, moreover, that φ ( u ) =Φ( | u | ) for some convex and nondecreasing C function Φ : R + → R + such that Φ(0) =0 = min w ∈ R + Φ( w ) . If the functions u and g in Lemma 4.1 are componentwise non-negative and g ∈ L (Ω; R m ) , then the minimizer u is (componentwise) nonnegative. Remark 5.1.
We can no longer allow g to be an H − (Ω; R m ) distribution. Physically,there cannot be any sinks, i.e. the sources must be non-negative functions f i = f i ( t, x ) ≥ for all t ≥ .Note also that the structural assumption on the potential φ can be weakened to φ ( u + ) − b ( u ) · u + ≤ φ ( u ) − b ( u ) · u , for fixed u ∈ R m + and for all u ∈ R m , which is a kind ofadditional convexity/monotonicity property. Proof.
It is enough to prove that, if u ≥ , then any minimizer u from Lemma 4.1remains (componentwise) non-negative. Indeed, iterating will give positivity u ( n ) ≥ at the discrete level for all n ≥ , thus the piecewise-constant interpolation u h ( t ) ≥ and the final solution u = lim h → u h ≥ as well.Our extra assumption on the potential readily implies that φ ( u + ) = Φ( | u + | ) ≤ Φ( | u | ) = φ ( u ) for all u ∈ R m , where u + denotes the componentwise positive part ( u + ) i = max( u i , (and of course | u + | ≤ | u | ). Moreover for all u ≥ we see that b ( u ) = Φ (cid:48) ( | u | ) u | u | ≥ hasnonnegative components as soon as u ≥ . As a consequence, it holds b ( u ) · u ≤ b ( u ) · u + . Then using u + in the minimization problem, we get E ( u + ) = (cid:90) Ω φ ( u + ) − b ( u ) · u + h dx + 12 (cid:90) Ω |∇ u + | dx − (cid:90) Ω g · u + dx ≤ (cid:90) Ω φ ( u ) − b ( u ) · uh + 12 (cid:90) Ω |∇ u | dx − (cid:90) Ω g · u dx = E ( u ) , where we have used Stampacchia’s result ∇ u + i = χ [ u i > ∇ u i and thus (cid:107)∇ u + i (cid:107) L ≤(cid:107)∇ u i (cid:107) L . Uniqueness of the minimizer implies that u + = u , from which it follows that u ( x ) ≥ componentwise and a.e. x ∈ Ω . (cid:3) Numerical results
In this section, we will complement the previous existence results by implementing asimple finite-difference algorithm to approximate the solution u of problem (4.1). Forthe ease of exposition we restrict ourselves to the one-dimensional case ( d = 1 ) with twospecies ( m = 2 ).Let us choose a fixed space interval [ a, b ] and, for fixed I ∈ N , we take a uniformpartition ∆ x = b − aI + 1 , x i = a + i ∆ x, i = 0 . . . I + 1 , so that x I +1 = b . Similarly, for a fixed T < ∞ and N ∈ N , we set ∆ t = TN , t n = n ∆ T, n = 0 . . . N.
The unknowns take values at the meshpoints ( x i , t n ) and, as usual, we write U ni ≈ u ( x i , t n ) , i = 1 , . . . , I , n = 0 , . . . , N . For simplicity, let us consider homogeneous Dirichlet boundary conditions u ( a, t ) = u ( b, t ) = 0 ⇔ U n = U nI +1 = 0 ∀ n , and assume that f = 0 . We construct a global approximation, continuous in x andpiecewise linear at each interval K i = ( x i , x i +1 ) , to approximate each of the terms inthe discrete system of equations. To this end, we first choose the elementary quadraturemethod (cid:90) x i +1 x i φ (cid:0) u ( x, t n +1 ) (cid:1) dx ≈ φ ( U n +1 i ) + φ ( U n +1 i +1 )2 ∆ x. Defining b ni = b ( u ni ) , we obtain similarly the approximation (cid:90) x i +1 x i b ( u ( x, t n )) · u ( x, t n +1 ) dx ≈ b ni · U n +1 i + b ni +1 · U n +1 i +1 x . Moreover, the Dirichlet energy is approximated as (cid:90) x i +1 x i |∇ u ( x, t n +1 ) | dx ≈ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) U n +1 i +1 − U n +1 i ∆ x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ∆ x. Following the approach presented in the previous sections, and after initialization U i = u ( x i ) , OMPETITIVE CONTAMINANT TRANSPORT 17 the approximate solution is determined by solving recursively(6.1) U n +1 = Argmin F n , where the real-valued function F n : ( R m ) I → R is the spatial approximation of F n ( u ) = (cid:90) Ω (cid:0) φ ( u ( x )) − b ( u n ( x )) · u ( x ) (cid:1) dx + ∆ t (cid:90) Ω |∇ u ( x ) | dx , and, according to the above-defined quadratures, is explicitly written as(6.2) F n ( U ) := I (cid:88) i =1 φ ( U i )∆ x − I (cid:88) i =1 b ni · U i ∆ x + ∆ t I (cid:88) i =0 (cid:12)(cid:12)(cid:12)(cid:12) U i +1 − U i ∆ x (cid:12)(cid:12)(cid:12)(cid:12) ∆ x . Note that both the left and right cells, [ x i − , x i ] and [ x i , x i +1 ] , contribute with weight / to the first two terms for each i . Moreover, we have exploited the zero boundaryconditions U = U I +1 = 0 in the last term and observed that no contributions φ (0) = 0 and b (0) = 0 arise from the first two terms for i = 0 , I + 1 . This construction is nothingbut the spatial discretization of the function u ( n ) defined (4.5) using the minimizationscheme of Lemma 4.1.All the simulations were performed using Gnu/Octave on a personal computer. Thenumerical minimization step can be performed using any unconstrained optimizationmethod. We have used the fminunc toolbox with gradient option, since the Jacobian D U F n can easily be computed analytically from the above formula. Remark 6.1.
We point out that our scheme could be adapted to arbitrary space di-mension using, e.g. , finite elements. One could also choose higher order quadratures,but the naive trapezoidal quadrature method used here has the advantage of allowing theJacobians to be computed analytically.
A two-component test-case.
As an example we present here a two-componentcomputation u = ( u , u ) for the Freundlich isotherm b ( u ) = u + | u | p − u , with exponent p = 1 / (for instance). The computation was performed for ∆ t = ∆ x = 10 − in adomain ( x, t ) ∈ ( − , × (0 , . (but other intervals could have been used) and theresult is shown in Fig. 1.It is worth pointing out that our scheme seems to be automatically positive, as shouldbe expected from the results of Section 5 since the particular isotherm b ( u ) = u + | u | / − u maps the positive cone ( R + ) into itself. In the same spirit and as alreadydiscussed, one should expect that the singularity/degeneracy D u b (0) ≈ ∞ leads tofinite speed of propagation and free boundary solutions. The scheme seems to capturereasonably well the propagation of free boundaries, as shown in Fig. 1; the left and rightfree-boundaries stay away from the boundaries x = ± for small times, and the internalhole at time t = 0 persists until t ≈ . (after which the supports of u ( t, . ) and u ( t, . ) match for all later times and progressively invade the whole domain with finite speed,which is a characteristic feature of degenerate diffusion).6.2. Error analysis.
In Section 4, we proved convergence of the semi-discretized so-lution u h to the unique weak solution u as the time-step h → , but without errorestimates. Therefore, any theoretical convergence result U → u with a priori error es-timates as ∆ x, ∆ t → for the fully discretized problem is certainly beyond the scopeof this paper. However, we wish to present instead a numerical investigation of theconvergence orders.To compare the numerical solution we need an analytical one. To the best of ourknowledge this is not possible in the two-component case. In the scalar case (singlecomponent) and for the specific isotherm b ( u ) = | u | p − u with Freundlich exponent p ∈ (0 , , the change of variables z = b ( u ) turns the homogeneous equation ∂ t b ( u ) = ∆ u into the celebrated Porous Medium Equation (PME, in short) ∂ t z = ∆( | z | m − z ) with m = 1 p > . As is well known [24], the Cauchy problem for PME is well posed in the whole space forsuitable initial data, and the problem exhibits finite speed of propagation, i.e. for anycompactly supported initial data z ( x ) the solution z ( t, . ) stays compactly supported for Figure 1. snapshot of u ( t, . ) (dashed line) and u ( t, . ) (solid line) asfunctions of x for several times from t = 0 to t = 0 . .all times, and the supports expand with finite speed. More importantly, the Zel’dovich-Kompaneets-Barenblatt profiles ( t, x ) (cid:55)→ Z ( C ; t ; x ; t, x ) := ( t + t ) − α (cid:0) C − k | x − x | ( t + t ) − β (cid:1) m − + ≥ yield a family of explicit solutions for t > − t and x ∈ R d . Here C > is a free parametercorresponding to the L ( R d ) mass (which is preserved along evolution), t ∈ R and x ∈ R d translate the invariance under shifts, and α = dd ( m −
1) + 2 , β = αd , k = α ( m − md are universal constants depending on the exponent m > and the space dimension d ≥ only. Moving back to u = | z | m − z yields explicit solutions ( t, x ) (cid:55)→ U ( C ; t ; x ; t, x ) tothe original equation, which we use to compute errors. The results from the previoussections, and in particular the a priori energy estimates, suggest that the implicit schemegiven by (6.1)-(6.2) should be stable unconditionally with respect to ∆ t , the stabilityoccurring in the relevant energy spaces corresponding to Theorem 1. For our tests wealways chose a fixed computation time T = . , and fix ∆ t = 2∆ x for convenience.The other parameters are adjusted so that the ZKB solutions stay supported inside thenumerical domain Ω = ( a, b ) , so that these free-boundary solutions correspond to theunique solution to the Cauchy problem in Ω with zero boundary conditions at least for t ≤ T .Note that the larger the diffusion exponent m > , the more degenerate the equation.The speed of propagation becomes infinite in the limit p = m ↓ since the PME ∂ t z = ∆ z m then formally converges to the heat equation, which has infinite speed ofpropagation as all uniformly parabolic equations. This is also captured by our scheme,but due to the lack of space we do not present the corresponding simulations here. The L ( Q T ) errors are shown in Fig. 2 as a function of ∆ x → for various values of m > ,suggesting algebraic convergence: the convergence seems to be affected only through aprefactor C m in e (cid:46) C m | ∆ x | r , but the rate r appears to be independent of m . OMPETITIVE CONTAMINANT TRANSPORT 19 −2.2 −2 −1.8 −1.6 −1.4 −1.2 −1−3−2.8−2.6−2.4−2.2−2−1.8−1.6−1.4 log ( ∆ x) l og ( e ) Error m=5m=2m=1.1
Figure 2.
Logarithmic plot of the L ( Q T ) errors e as a function of ∆ x , with ∆ t = 2∆ x . Acknowledgments.
This work was partially supported by FCT/Portugal through theproject UID/MAT/04459/2013 and by the UT Austin|Portugal CoLab project. FB andLM were supported by the Portuguese National Science Foundation through throughFCT fellowships SFRH/BPD/ 33962/2009 and SFRH/BPD/88207/2012.
References [1] Aizinger V., C. Dawson, B. Cockburn and P. Castilho, The local discontinuous Galerkin method forcontaminant transport,
Adv. Water Res. (2001), 72-78.[2] Alt W. H and S. Luckhaus, Quasilinear Elliptic-Parabolic Differential Equations, Math. Z. (1983),311-341.[3] Amundson N.R, R. Aris and H.K Rhee,
First-order Partial Differential Equations , Englewood Cliffs,N.J, Prentice-Hall, 1989.[4] Barrett, J.W. and P. Knabner, Finite element approximation of transport of reactive solutes in porousmedia. Part 1. Error estimates for non-equilibrium adsorption processes,
SIAM J. Numer. Anal. (1997), 201-227.[5] Barrett J. W. and P. Knabner, Finite element approximation of the transport of reactive solutes in porousmedia. Part II: Error estimates for equilibrium adsorption processes, SIAM J. Numer. Anal. (1997),455-479.[6] Bear J., Dynamics of Fluids in Porous Media , Elsevier, New York, 1972.[7] Bear J. and A. H-D. Cheng,
Modeling Groundwater Flow and Contaminant Transport , Springer, NewYork, 2010.[8] Daskalopoulos, P and C.E. Kenig,
Degenerate Diffusions: Initial Value Problems and Local RegularityTheory , EMS, Zürich, 2007.[9] Dawson, C., van Duijn, C.J. and R.E. Grundy, Large-time asymptotics in contaminant transport inporous media,
SIAM J. Appl. Math. (1996), 965-993.[10] Dawson, C., van Duijn, C.J. and M.F. Wheeler, Characteristic-Galerkin methods for contaminant trans-port with non-equilibrium adsorption kinetics, SIAM J. Numer. Anal. (1994), 982-999.[11] Fonseca I. and G. Leoni, Modern Methods in the Calculus of Variations with Applications to NonlinearContinuum Physics , Springer, 2007.[12] Gutierrez, M. and H.R. Fuentes, Modeling adsorption in multicomponent systems using a Freundlich-typeisotherm,
J. Contam. Hydrol. (1993) 247-260.[13] Hinz, C., Description of sorption data with isotherm equations, Geoderma (2001), 225-243.[14] Kuusi, T., Monsaingeon, L. and J.H. Videman, Systems of partial differential equations in porous medium, Nonlinear Anal. (2016), 79-101.[15] Limousin, G., Gaudet, J.-P., Charlet, L., Szenknect, S., Barthès, V. and M. Krimissa, Sorption isotherms:A review on physical bases, modeling and measurement,
Appl. Geochem. (2007), 249-275.[16] Monsaingeon, L., An explicit finite-difference scheme for one-dimensional Generalized PorousMedium Equations: interface tracking and the hole filling problem, ESAIM: M2AN , in press.http://dx.doi.org/10.1051/m2an/20150630 BAÍA, BOZORGNIA, MONSAINGEON, AND VIDEMAN[17] Rockafellar R. T.,
Convex Analysis , Princeton University Press, 1997.[18] Rothe, E., Zweidimensionale parabolische Randwertaufgaben als Grenzfall eindimensionaler Randwer-taufgaben,
Math. Ann. (1930), 650-670.[19] Roubíček T.,
Nonlinear Partial Differential Equations with Applications , Birkhäuser Verlag, 2005.[20] Sheindorf, C., Rebhum, M. and M. Sheintuch, A Freundlich type multicomponent isotherm,
J. Colloid.Interface Sci. (1981), 136-142.[21] Srivastava, V.C., Mall, I.D. and I.M. Mishra, Equilibrium modelling of single and binary adsorption ofcadmium and nickel onto bagasse fly ash, Chem. Eng. J. (2006) 79-91.[22] Susarla, S., Bhaskar, G.V. and S.M.R. Bhamidimarri, Competitive adsorption of phenoxy herbicide chem-icals in soil, Water Sci. Tech.
26 (1992), 2121-2124.[23] Wu C-H., C-Y. Kuo, C-F. Lin and S-L. Lo, Modeling competitive adsorption of molybdate, sulfate,selenate and selenite using a Freundlich-type multi-component isotherm,
Chemosphere (2002), 283-292.[24] Vázquez, J.L., The Porous Medium Equation: Mathematical Theory , Oxford University Press, NewYork, 2007.
CAMGSD/Departamento de Matemática, Instituto Superior Técnico, Universidade deLisboa, Av. Rovisco Pais 1, 1049-001 Lisboa, Portugal
Email address : [email protected] CAMGSD, Instituto Superior Técnico, Universidade de Lisboa, Av. Rovisco Pais 1,1049-001 Lisboa, Portugal
Email address : [email protected] CAMGSD, Instituto Superior Técnico, Universidade de Lisboa, Av. Rovisco Pais 1,1049-001 Lisboa, Portugal
Email address : [email protected] CAMGSD/Departamento de Matemática, Instituto Superior Técnico, Universidade deLisboa, Av. Rovisco Pais 1, 1049-001 Lisboa, Portugal
Email address ::