Variational – hemivariational system for contaminant convection – reaction – di ﬀ usion model of recovered fracturing ﬂ uid

: This work is devoted to study the convection – reaction – di ﬀ usion behavior of contaminant in the recovered fracturing ﬂ uid which ﬂ ows in the wellbore from shale gas reservoir. First, we apply various constitutive laws for generalized non-Newtonian ﬂ uids, di ﬀ usion principles, and friction relations to formulate the recovered fracturing ﬂ uid model. The latter is a partial di ﬀ erential system composed of a nonlinear and nonsmooth stationary incompressible Navier-Stokes equation with a multivalued friction boundary condition, and a nonlinear convection – reaction – di ﬀ usion equation with mixed Neumann boundary conditions. Then, we provide the weak formulation of the ﬂ uid model which is a hemivariational inequality driven by a nonlinear variational equation. We establish existence of solutions to the recovered fracturing ﬂ uid model via a surjectivity theorem for multivalued operators combined with an alternative iterative method and elements of nonsmooth analysis.


Introduction
The study is motivated by recent models in the production of the shale gas.Shale gas is a natural gas that is located in some geological formations and is considered today as an important unconventional natural gas resource.The US, Canada, and P.R.China are the only few countries producing shale gas for commercial use.To overcome an extremely low permeability of shale, for the production process, hydraulic fracturing techniques have been developed.They consist of pumping fracturing fluid into the wellbore under a high pressure in order to form artificial fractures and create a natural gas outflow channel.A mixture of the fracturing fluid and the formation water is usually called the recovered fracturing fluid (or the recovered fluid).Before placing a well on production, the fracturing fluid must be pumped out of the wellbore (in the flowback process) since the contaminants of recovered fracturing fluid may affect the surface soil, surface water, and shallow groundwater.During the flowback process, the data can be gained on hydraulic fractures, the effectiveness of the fracturing operations, as well as on the prediction of long-term reservoir performance.The studies on the recovered fracturing fluid and the diffusion of contaminants in the hydraulic fracturing processes attracts attention of many researchers.We refer, for instance, to [4,10,21,30,35] and references therein for more details.
In this work, we study a system of a hemivariational inequality for the stationary incompressible Navier-Stokes equation with a multivalued friction boundary condition and a nonlinear convection-reaction-diffusion equation with mixed Neumann boundary conditions.The system is coupled by a nonconvex and nondifferentiable superpotential which appears in the nonmonotone slip boundary condition, see (3.11), and which depends on the concentration of the contaminant in recovered fracturing fluid.The second coupling member is the convection term, see (3.18), which appears in a nonlinear convection-reaction-diffusion equation and depends on the velocity field of the recovered fracturing fluid.To the best of our knowledge, this formulation establishes a new mathematical model to study simultaneously the flow behavior of the recovered fracturing fluid and the reaction-diffusion phenomenon of contaminants in the wellbore of shale gas reservoir.
The main result of the study is an existence theorem for weak solutions to the aforementioned system, see Theorem 5.1.The method of proof is based on a surjectivity theorem for multivalued pseudomonotone operators combined with an alternative iterative approach.The theorem contains a smallness hypothesis on the constants.Removal of such condition is an interesting open problem.Moreover, we mention that it is also of interest to study an evolution counterpart of the system.
Note that the notion of a hemivariational inequality is based on the Clarke generalized subgradient introduced for locally Lipschitz functionals.In comparison with the variational inequalities which emerge from convex energy principles, the stationary hemivariational inequalities lead to substationarity problems and are closely related to nonsmooth and nonconvex energies.The investigation on the hemivariational inequalities has been originated in 1980s with the works of Panagiotopoulos [31,32].For the recent developments, we refer to some representative results [1,3,11,[13][14][15]18,20], monographs [6,7,26,29,34], and references cited therein.The mathematical models of the steady state flow of the incompressible fluid with various types of multivalued and nonmonotone boundary conditions have been studied only recently in [8,9,16,19,[22][23][24][25]27,36].The analysis of these fluid flow problems leads in a natural way to variational or variational-hemivariational inequalities for the velocity field.
The outline of this study is as follows.Section 2 collects the necessary notations and preliminary results.In Section 3, we establish the recovered fracturing fluid model to describe the convection-reaction-diffusion behavior of contaminant in the recovered fracturing fluid which flows in the wellbore from shale gas reservoir.Then, in Section 4, we impose the mild assumptions of the recovered fracturing fluid and obtain its weak formulation, which is a variational-hemivariational system, see Problem 4.2.Finally, in Section 5, an existence result to the recovered fracturing fluid model under consideration is established.

Mathematical prerequisites
In this section, we shortly review a preliminary material from convex and nonsmooth analysis and recall the basic notations. Let The following result collects properties of the generalized directional derivative and the generalized subgradient in the sense of Clarke (see, for instance, [26,Proposition 3.23]), which is used in the subsequent sections.Let V be a reflexive Banach space with its dual space V * and → F V : 2 V * be a set-valued map.Then, F is called to be pseudomonotone in the sense of Brézis, if the following conditions are satisfied (i) for every ∈ v V, ( ) F v is nonempty, bounded, closed, and convex in V *, (ii) F is u.s.c.from every finite-dimensional subspace of V to V * endowed with the weak topology, We conclude this section by recalling the following surjectivity theorem for multivalued mappings, which is formulated for the sum of a maximal monotone multivalued operator and a bounded multivalued pseudomonotone mapping.The proof of this theorem can be found in [17,Theorem 2.2].
Theorem 2.2.Let V be a real reflexive Banach space, 2 V * be a bounded set-valued pseudomonotone operator, and ∈ f V*.Assume that there are 3 Recovered fracturing fluid model end, we apply various constitutive laws for generalized non-Newtonian fluids, diffusion principles, and friction relations to obtain its classical formulation.We assume that the recovered fracturing fluid model is considered in the following physical setting.The horizontal wells have no influence on each other, which means that the formation of a single horizontal well and a pressure drop between the wellbore and the shale gas reservoir will not be influenced by each other (Figure 1).Additionally, the recovered fracturing fluid in the model is supposed to be non-Newtonian, the interior of the horizontal wellbore does not have porous medium, and the flow state of recovered fracturing fluid would be regarded to be the free flow.Therefore, from the physical point of view, the flow behavior of recovered fracturing fluid fulfills an incompressible Navier-Stokes equation for non-Newtonian fluid.Moreover, we suppose that • the chemical contaminants are diffused in the recovered fracturing fluid, which is a non-Newtonian fluid, so, anomalous diffusion happens in our model (namely, generalized nonlinear Fick's diffusion law is satisfied); • the velocity field to recovered fracturing fluid causes the presence of the convection term ⋅∇ v c which involves the velocity field v of the recovered fracturing fluid and the gradient ∇c of the concentration of the contaminant; • the chemical contaminants in the recovered fracturing fluid may lead to a chemical reaction with the minerals and a catalytic agent in horizontal wellbore (thus, a reaction term ( ) x g c , will be involved in the model).
Therefore, the concentration of the contaminant in the recovered fracturing fluid enjoys a nonlinear convection-reaction-diffusion equation.To conclude, we can see that the recovered fracturing fluid model under consideration is described by a partial differential system that couples a nonlinear Navier-Stokes equation and a nonlinear convection-reaction-diffusion equation.2).Throughout the study, we often do not indicate explicitly the dependence of various functions on the spatial variable ∈ = ∪ x Ω Ω Γ.We use the symbols described in Table 1.Next, we move our attention to describe the recovered fracturing fluid model.The mass density ρ of the recovered fracturing fluid is supposed to be well-defined.For any given subregion ⊂ B Ω, we suppose that it has smooth boundary ∂B.Let n be the unit outward normal defined on the boundary ∂B and v be the velocity field of the recovered fracturing fluid in Ω.It is not hard to see that ⋅ v n and ⋅ v ρ n represent the volume flow rate per unit area on the boundary ∂B and the mass flow rate per unit area, respectively.The law of conservation of mass reveals that the identity is true: In other words, the following integro-differential equation holds Directly by the divergence theorem, we have Hence, we obtain the following integral formulation of the law of conservation of mass B Ω is arbitrary, we use the above identity to obtain the following continuity equation: Moreover, we assume that the rate of change of mass density is quite small with respect to time and space variables, thus, ρ is a constant.For simplicity, in what follows, we take ≡ ρ 1.Therefore, (3.1) can be rewritten as the following identity: The inner product in d , i.e., = ξ τ ξ τ : The latter is usually called the divergence-free condition, which means that the recovered fracturing fluid is incompressible.Let ( ) = = x x i i d 1 be the path followed by a recovered fracturing fluid particle.It is well-known that the velocity field v and acceleration a of the fracturing fluid particle fulfill the following equalities: Then, we have where ⋅∇ v v is given by Next by π we denote the pressure of the recovered fracturing fluid.Let ( ) ( ) , 1 be the deformation tensor of velocity field v with ( ) for all , 1,…, .
If the recovered fracturing fluid is a Newtonian fluid, it follows from the fact that the stress tensor of a moving fluid approaches the stress tensor of a stationary fluid when the motion stops.This means that the stress tensor P could be formulated as the sum of the isotropic part − I π and the anisotropic part where the last inequality is obtained by using the divergence-free condition (3.2).Here ∈ I d is the identity matrix, μ is the first coefficient of viscosity, and ζ is the second coefficient of viscosity.Let × → C : Ω d d be a Cauchy stress tensor which is, in general, nonlinear, because the recovered fracturing fluid is considered to be a non-Newtonian fluid.Therefore, it is reasonable to assume that the stress tensor P satisfies the following generalized constitutive law: where ( ( )) , the fluid is Newtonian), then (3.5) reduces to the classical law (3.4).In what follows, on the boundary Γ, we denote by , the traction vector of the total stress tensor.So, are the normal and tangential components of the traction vector ( ) τ v π , on the boundary, respectively.It is clear that Keeping in mind the fact that the force of stress exerted across the surface ∂B per unit area is Hence, for any vector ∈ ξ d fixed, we use the divergence theorem to find The arbitrariness of ∈ ξ d reveals that We denote by f the given body force applied per unit mass of the recovered fracturing fluid in Ω.Then, the total body force acted in B is expressed by

B
Hence, the total force of the recovered fracturing fluid in the subdomain B is given by Recall that ⊂ B Ω is arbitrary, thereby, we have By the Newton second law, we obtain the following equation of balance of momentum: Inserting (3.3) and (3.6) into the equality above and using ≡ ρ 1, it yields the following nonlinear Navier-Stokes equation: In this study, we are interested in the study of the stationary counterpart of the nonlinear Navier-Stokes equation (3.7) of the form: Next we describe the boundary conditions for the recovered fracturing fluid.Since the wellbore is considered to be half-closed, so, from the mathematical point of view, it is reasonable to suppose that the recovered fracturing fluid is adhered to the wall, say, the boundary Γ D .This means that the velocity field v of the recovered fracturing fluid satisfies the homogeneous Dirichlet condition on Γ D , i.e., = v 0 on Γ .

D (3.9)
Let Γ C be the boundary such that there is no phenomenon of osmosis to the recovered fracturing fluid, that is, there is no inflow and outflow on the boundary Γ C .Therefore, it is not difficult to see that the normal component of velocity field to the recovered fracturing fluid satisfies the no-flux boundary condition on Γ C , i.e.
On the other hand, the tangential component of the traction vector ( ) τ v τ and the velocity field v are assumed to obey the following multivalued and nonmonotone friction law on boundary Γ C , namely, , , , on Γ , where c stands for the concentration of the contaminant in the recovered fracturing fluid, and , , is locally Lipschitz continuous for a.e.∈ x Γ C and all ∈ s , and ( ) ∂ x ξ j s , , is the generalized Clarke subgradient of j with respect to its last variable.Condition (3.10) is called the impermeability (no leak) boundary condition, and condition (3.11) is called the nonmonotone slip boundary condition.In particular, if j is of the form C is a given function, then the multivalued friction condition (3.11) can be written as follows: By the definition of the convex subgradient, we know that the above inclusion can be rewritten equivalently as the following Tresca friction law (see [26,Section 6.3] for a detailed discussion): It should be pointed out that in this study, we suppose that the potential j is nonconvex and locally Lipschitz with respect to its last variable.Furthermore, we assume that the domain Γ I connects the wellbore and fractures such that the recovered fracturing fluid and contaminants can flow in/out from fractures to the wellbore.For simplicity, we only consider the situation that the recovered fracturing fluid flows in from fractures to wellbore.On boundary Γ I , we suppose that the following generalized Signorini-type contact condition holds: where ≥ ϱ 0 is a constant and → ϕ : Γ I is the measured force on the boundary Γ I .Such boundary condition means that the recovered fracturing fluid and contaminants are only injected into wellbore from fractures.Also, in order to prolong the service life of wellbore, a restraint mechanism to the normal component of velocity field v ν is considered, i.e., the rate of inflow of the velocity field cannot be able to exceed ϱ (the negativity stands for inflow effect) and the Signorini-type contact condition for the normal velocity (i.e., complementary condition (3.12) 2,3 ) is considered.Assume that there is no friction in the tangential direction for the recovered fracturing fluid on the outflow boundary Γ O .So, from the incompressibility of the recovered fracturing fluid, it is not difficult to see that the recovered fracturing fluid satisfies the following constitutive laws on the outflow boundary where → φ: Γ O is a measurable function which is generated by the surface force.In the sequel, we denote by c the concentration of contaminants in the wellbore.If the diffusion behavior of contaminants satisfies the well-known Fick diffusion law, that is, the rate of diffusion of a substance across unit area is proportional to the concentration gradient (see, for instance, [33]), then it holds where E is the molecular diffusion flux of contaminants and κ is a diffusion coefficient.When the coefficient of diffusion of contaminants is affected by the velocity field of the recovered fracturing fluid, then the Fick diffusion law (3.14)can be generalized to the following one: . Since the recovered fracturing fluid is considered to be non-Newtonian and anisotropic, the diffusion of contaminants satisfies the non-Fickian diffusion law, namely, the following nonlinear and nonhomogeneous diffusion equation holds Variational-hemivariational system for contaminant convection-reaction-diffusion model  9 The latter combined with the law of mass conservation implies Moreover, if some chemical reactions take place in the wellbore, then a certain reaction term × → g : Ω will appear in the diffusion equation above, i.e., the reaction-diffusion equation reads as follows: Furthermore, when the velocity field v influences the diffusion of contaminants (i.e., the presence of convec- tion), then the above reaction-diffusion equation can be modified to the following convection-reaction-diffusion equation: (3.17) As we have mentioned before, we will focus on the investigation of the stationary diffusion process of the contaminants and the flow of the recovered fracturing fluid in the wellbore.We arrive at the convectionreaction-diffusion equation Recall that the contaminants could flow in from the fractures to the wellbore only on the boundary Γ I and flow out from wellbore to outside on the boundary Γ O , so the following mixed Neumann boundary conditions are assumed: • there is no flux of contaminants through the boundary ∪ Γ Γ D C , thus, • the contaminants enter (or expel) from the wellbore through fractures on the boundary Γ I (or Γ O ) which means that is a given source term on the inflow-outflow boundary ∪ Γ Γ

Hypotheses and the weak formulation
The goal of this section is two-fold.The first one is to deliver mild hypotheses on the data of the recovered fracturing fluid model, Problem 3.1, which will be used in Section 5 to prove the existence of weak solutions to the recovered fracturing fluid model.The second goal is to apply the variational analysis technique and the properties of the Clarke subgradient to establish the weak formulation of Problem 3.1, which is a coupled system consisting of a highly nonlinear hemivariational inequality with constraints and a nonlinear variational equation.First, we impose the following hypotheses on the data of Problem 3.1.

( )
is continuously differentiable (i.e., C 1 ) and strictly convex for a.e.∈ x Ω with ( ) (iii) there exists a function for all ∈ D d and a.e.∈ x Ω; (iv) the inequality Variational-hemivariational system for contaminant convection-reaction-diffusion model  11 holds for all ∈ D d and a.e.∈ x Ω with ( ) ( ) , , is locally Lipschitz for all ∈ s and a.e.∈ x Γ C ; (iii) there exist a function for all , and a.e.Γ , where ∂j stands for the generalized Clarke subgradient operator of j with respect to its last variable; , , ; ( ) x Ω; (iii) there exist two constants > a b , 0 , for all and a.e.Ω.
x Ω; (iii) there exist a function where p * 2 represents the Sobolev critical exponent of p 2 ; (iv) there exist a function ( ) , and ≥ ϱ 0. There are plenty of functions which satisfy ( ) H g .A concrete example is as follows ( To give a weak formulation of Problem 3.1, we introduce the subspace of ( ) We denote by E, the closure of in ( ) W Ω; Since, the measure of Γ D is positive, it permits us to invoke Korn's inequality (see, for instance, [25]) to find a constant where ( ) D v is the deformation tensor of v. Hence, we note that Applying the Green formula to the second term on the left-hand side of (4.1), we where we have used the boundary conditions = = w v 0 on Γ D , = = w v 0 ν ν on Γ C , and the divergence-free condition for w and v. Also, we use the divergence theorem (see, e.g., [26,Theorem 2.25]) to the differential operator By the definition of ( ) τ v π , and (3.5), we have Recall that so we use the boundary condition (3.24) to obtain Variational-hemivariational system for contaminant convection-reaction-diffusion model  13 Therefore, we obtain Applying the boundary conditions (3.25)-(3.27)and the definition of the generalized subgradient in the sense of Clarke, it yields for all ∈ w K .By an analogous argument applied to the convection-reaction-diffusion equation (3.28), we derive the following variational equation: Taking into account the above equation and (4.3), we are now in a position to give the weak formulation of Problem 3.1 as follows.

Main existence result
The main goal of this section is to establish an existence theorem for weak solutions to Problem 3.1.The main method is based on a surjectivity theorem for pseudomonotone operators together with an alternative iterative approach.The existence theorem for weak solutions to Problem 3.1 is stated as follows.
, where the function ( ) { } +∞ → δ : 0, 0, 1 is defined by In order to prove this theorem, we need the following result.Consider the nonlinear operator for v, ∈ w E. The following lemma provides several critical properties of the operator A, which will be applied to establish the existence of solutions to Problem 4.2.

Lemma 5.2. If hypothesis ( ) C H
holds, then the operator → A E E : * defined in (5.1) is well-defined, bounded, continuous, maximal monotone, and of type ( + S ).
Proof.Let ∈ v E be arbitrary.Using the hypothesis ( ) C H (iii), Hölder's inequality, and the elementary inequality ( ) ( ) for some > C 0 0 .This means that A is well-defined, bounded, and satisfies the growth condition is continuously differentiable and strictly convex for a.e.∈ x Ω, so it is easy to see that Then, by the Lebesgue dominated conver- gence theorem, we have This proves that A is continuous on E. Employing [2, Theorem 1.3, p. 45], we conclude that A is maximal monotone.It remains to show that A is of ( Hence, we have The monotonicity of A points out and hence . Passing to a subsequence, if necessary, we may assume that n By the convexity of G, and the fact ( for a.e.∈ x Ω and all ∈ D d .The assumption H(C)(iv) guarantees that for all ∈ D d and a.e.∈ x Ω.The latter combined with hypothesis ( ) Therefore, we are able to find a function with where It follows from hypotheses ( ) H j , [26, Theorem 3.47], and [28, Lemma 7] that functional defined by (5.6) satisfies the following properties: , , ; dΓ . Moreover, when = d 3 and < p 3 1 , the embedding of E into ( ) is compact and Therefore, we use the definition of B, the generalized Hölder inequality, and the Korn inequality to find a constant > c 0 B such that (see, for example, [12, Lemma 2.1, p. 284]) The latter combined with the definition of the form b in (5.5) implies that B is linear and continuous with respect to each variable.Next we shall show that [ ] .
Passing to the limit as → ∞ n in the estimate above, we have c Variational-hemivariational system for contaminant convection-reaction-diffusion model  17 We assert that c is pseudomonotone.Based on the above analysis, it is obvious that ( ) ⋅ c is bounded and for each ∈ v E, the set ( ) v c is nonempty, closed, and convex in E*.So, from [26, Proposition 3.58], it is sufficient to prove that c is generalized pseudomonotone.
Let Hence, for each ∈ n , there is an element 2 for a.e.∈ x Γ.Next for each ∈ w E, we use the Green formula and the divergence-free condition to derive Letting → ∞ n in the above equality and using the Lebesgue dominated convergence theorem, we obtain The latter combined with the weak-weak continuity of B implies Recall that A has the ( + S )-property (see Lemma 5.2), so, we infer that and the continuity of the embedding of E into ( ) L Γ; Passing to a subsequence, if necessary, we may suppose that there is ∈ ξ X* such that → ξ ξ n w in X*.Further, we apply the strong-weak upper semicontinuity of ( ) . This means that c is a generalized pseudomonotone operator.So, applying [26, Proposition 3.58], we conclude that c is pseudomonotone.
Furthermore, we shall show that c is coercive.For any ∈ v E fixed, by ( ) C H (ii) and (iv), one has The boundary conditions (3.26) and (3.27) show that By the definition of the generalized Clarke subgradient and hypothesis ( ) H j (iii), we have Taking into account the estimates above, we obtain Next we use the inequality that is, c is coercive.Using a standard procedure (for example, the Hahn-Banach theorem), we can show that ∈ v K is a solution to the following hemivariational inequality if and only if it solves the inclusion problem where I K is the indicator function of K , ∂ I C K is the convex subdifferential operator of I K , and ͠ ∈ f E* is defined by which is obtained via using Riesz's representation theorem.Since ∈ K 0 , one has where the constants ≥ c d , 0 defined by , is continuous and the following inequality holds Using the boundary conditions (3.24)-(3.27)and the divergence-free condition for the velocity field, we obtain where γ 2 is the trace operator from . We use conditions ( ) H κ and ( ) Because of ≥ θ p 2 and the smallness condition . Then, we have Subsequently, we apply the embedding theorem, the Lebesgue dominated convergence theorem, and hypothesis ( ) H g to obtain Hence, we have So, by the ( + S )-property of the operator . This means that H is pseudomonotone.We are now in a position to invoke [26, Theorem 3.74] and are both uniformly bounded.We only show that the sequence { } v n is uniformly bounded in E, since the boundedness of { } ( ) can be obtained directly by using a similar argument.Suppose that { } v n is unbounded.Without loss of generality, we may assume that ‖ ‖ → +∞ v n E as → ∞ n .For any ≥ n 1, we insert = w 0 into (5.9) to find we obtain a contradiction  .Finally, we pass to the upper limit as → ∞ n in (5.9) and (5. .Consequently, we conclude that ( ) v c , is a solution to Problem 4.2.This completes the proof of the theorem.□

Figure 1 :
Figure 1: The structure diagram of a wellbore.

3 .
Let ⊂ Ω d be a bounded domain (open and connected set) occupied by the viscous incompressible and non-Newtonian fluid (recovered fracturing fluid) with = d 2 and = d The boundary = ∂ Γ Ω is Lipschitz and Γ is partitioned into four mutually disjoint parts Γ C , Γ D , Γ I , and Γ O with Γ D having positive measure (Figure

dd
is continuous and monotone for a.e.∈ x Ω.Hence, it holds

C
By virtue of hypotheses ( )C H (ii), (iv) and ( ) H j (iii), and the boundary condition (3.26), we havePassing to the upper limit as → ∞ n in the inequality above, and using the inequalities ≥

C( 2 ( 5 . 14 )where we have used the monotonicity of the p 2 -
upper limit as → ∞ n in (5.11), and using (5.12) and (5.13), we obtain The latter together with the weak convergence → v v n w in E and Lemma 5.2 implies that → v v n in E. Next putting = − z c c n in (5.10), we have It follows from the Sobolev embedding theorems, hypothesis ( ) H g , and the Lebesgue dominated convergence theorem thatUsing the latter and taking the upper limit in(5.14)  as → ∞ n Laplace operator and hypothesis ( ) H κ (iii).Therefore, we apply the ( + S )-property of the operator ( > p 1 and ≥ d 2. In the sequel, we denote by ′ p and p*, the conjugate exponent and Sobolev critical exponent of p, respectively, i.e., "), the dual space of E, the duality pairing of E* and E, and the strong (weak) convergence in E, respectively.We say that a function E , and "→" (resp."→ w

Table 1 :
Symbol description otherwise .Taking into account (3.2), (3.8)-(3.13),(3.18), and (3.21), we have the classical formulation of the recovered fracturing fluid model.It is described by a nonlinear and nonsmooth stationary incompressible Navier-Stokes equation with a multivalued friction boundary condition coupled with a nonlinear convection-reaction-diffusion equation with mixed Neumann boundary conditions.
in what follows, we adopt ‖ ‖ ⋅ E as the norm of E.Moreover, it is easy to prove that E equipped with the norm ‖ ‖ ⋅ E becomes a separable and reflexive Banach space.In what follows, we denote by Proof of Theorem 5.1.The proof is divided into four steps.
n Furthermore, the boundedness of { ( ) Because of the maximal monotonicity of ∂ IC K , we observe that all conditions of Theorem 2.2 are verified.Using this theorem, we conclude that there exists ∈ v K which solves the inclusion problem(5.8).By the definition of the convex subdifferential and the generalized Clarke subgradient, we can see that ∈ c to obtain that, for each ∈ v K , equation (4.5) has at least one solution.
Variational-hemivariational system for contaminant convection-reaction-diffusion model  21The weak-weak continuity of B implies C(5.11)