A nonlinear model for marble sulphation including surface rugosity: theoretical and numerical results

We consider an evolution system describing the phenomenon of marble sulphation of a monument, accounting of the surface rugosity. We first prove a local in time well posedness result. Then, stronger assumptions on the data allow us to establish the existence of a global in time solution. Finally, we perform some numerical simulations that illustrate the main feature of the proposed model.


1.
Introduction. Deterioration and damage of stones is a complex problem and one of the main concern for people working in the field of conservation and restoration of cultural heritage. It is extremely difficult to isolate a single factor in this kind of processes, which are the results of the interaction of various mechanisms. These processes can be described through free boundary models, as in the case of damage induced by pollution, or by using a phase field approach. We introduce a new differential system coupling bulk and surface evolution equations to describe the phenomenon of marble sulphation of a monument, including the surface rugosity.
More precisely, we consider a monument made of calcium carbonate stone, located in a bounded domain Ω Ă R 3 , with a smooth boundary Γ, and subjected to a degradation process during a time interval p0, T q, for any given T ą 0. Following the approach introduced in [1,5,8,9,10,11] for related models of marble sulphation, we consider the system pφpcqsq t´d iv pφpcq∇sq "´λφpcqsc, in Ωˆp0, T q, (1.1) c t "´λφpcqsc, in Ωˆp0, T q, (1.2) φpcqB n s "´νprqps´sq, on Γˆp0, T q, (1.3) r t`B W prq`Ψ 1 prq`Gpr, c, sq Q F, on Γˆp0, T q, (1.4) equipped with the set of initial conditions sp0q " s 0 , in Ω, cp0q " c 0 in Ω, rp0q " r 0 on Γ. (1.5) Here s stands for the SO 2 porous concentration inside the material, c represents the local density of CaCO 3 and r ě 0 denotes the rugosity of the surface. We indicate by n the outward normal derivative to Γ. The first equation describes the evolution of the porous concentration of SO 2 , where φpcq is the porosity of the medium. This evolution is driven by the Fick's law and on the right hand side we have the reaction term between SO 2 and calcium carbonate, with rate λ ą 0. The second equation takes into account for the loss of calcium carbonate due the reaction.
The novelty in this model with respect to [1,8] is the introduction of the Robin boundary condition (1.3) for the flux of SO 2 through the external boundary, depending not only on the porosity of the medium and the external concentration s, but also on an effective permeability coefficient ν, which is a function of the superficial rugosity r.
Let us now describe what r represents. From the physical point of view the rugosity is quite a complex quantity, corresponding locally to the microscopic variation of a surface with respect to a flat configuration. For a precise definition see, for instance, [16,17]. Actually, in our model r stands for an effective macroscopic parameter which is formally a local damage parameter. As r grows, the profile of the surface has a greater deviation from the flat configuration and the microscopic density of peaks and valleys increases. In practice, this corresponds to a greater surface exposed to the pollution action, implying an increasing of the SO 2 flux permeability coefficient. Here, the value r " 0 corresponds to a completely smooth surface. Moreover, we deal with a material in which the "direction" of the evolution of r is not a priori fixed, i.e. r may increase or decrease through different external actions or repairing itself. The evolution of r is governed by the damage-like equation (1.4) (cf., e.g., [4,7,13] and the references therein). The symbol BW stands for the subdifferential of a proper, convex, and lower semicontinuous function W , accounting for possible internal constraints on r. The function G depends in particular on c and s through the porosity function φpcq and satisfies Gpr, 0, sq " Gpr, c, 0q " 0. The function Ψ 1 is sufficiently smooth and accounts for the non-monotone dynamics of r. The function F on the right hand side of (1.4) denotes possible external actions responsible for the formation of rugosity on the boundary, like, e.g., wind, rain or temperature variations.
It is worth noting that, as many dissipative evolution equations, our system (1.1)-(1.4) has a gradient flow structure (see, e.g., [12] where p G is such that G " B r p G. Then, using the kinetic equation (1.2) and letting Actually here we do not exploit this structure since we use only standard energy arguments, but this remark could be helpful for future developments.
In this paper we first prove a local in time well posedness result for system (1.1)-(1.5) in finite energy spaces. Next, under slightly stronger assumptions, we can establish some uniform bounds for the solution which in turn imply the global in time existence. Finally, we perform some numerical finite element simulations to assess the behavior of our model according to different physical situations.
The plan of the paper is the following. Section 2 is devoted to the statement of the main theorems and assumptions. The proof of the local existence result is detailed in Section 3. Section 4 contains the proof of the global existence result. Finally, in Section 5 we present the numerical simulations.
2. Main results. In this section, we introduce the formulation of system (1.1)-(1.5) in the correct functional setting and we state the main results. From now on, for the sake of simplicity, we will take λ " 1.
First of all, we assume that the porosity function φ is a linear function of the density c, i.e., φpcq " A`Bc, (2.1) where A and B are constants (see [9] and the references therein). Concerning the data, throughout the paper we make use of the following assumptions W : r0,`8q Ñ r0,`8s is convex and l.s.c., W p0q " 0, r 0 P L 2 pΓq, W pr 0 q P L 1 pΓq, r 0 ě 0, a.e. in Γ, (2.7) φpc 0 qB n s 0 "´νpr 0 qps 0´s q, a.e. in Γ, (2.8) We recall that if W is convex and lower semicontinuous then the subdifferential BW is a maximal monotone operator defined as BW pσ 0 q " tξ P R : W pσq ě W pσ 0 q`ξpσ´σ 0 qu.

Remark 2.3.
It is reasonable to impose a constraint on r so that r P r0, R 0 s, where R 0 ą 0 depends on the crystals dimension for the material we are considering. Actually, we can ensure the validity of this internal constraint assuming in (2.2), e.g., We can now formulate our problem Problem pP q. Find a triplet ps, c, rq such that and satisfying pφpcqsq t´d iv pφpcq∇sq "´φpcqsc, a.e. in Ωˆp0, T q, (2.15) φpcqB n s "´νprqps´sq, a.e. on Γˆp0, T q, (2.16) c t "´φpcqcs, a.e. in Ωˆp0, T q, (2.17) r t`ξ`Ψ 1 prq`Gpr, c, sq " F, ξ P BW prq, a.e. on Γˆp0, T q, (2.18) sp0q " s 0 , cp0q " c 0 , a.e. in Ω, rp0q " r 0 , a.e. on Γ. (2.19) We can now state our main existence results. The first theorem concerns the existence of a local in time solution.
Theorem 2.4. Let T ą 0 and (2.2)-(2.9) hold. Then there exists a time p T P p0, T s such that Problem pP q admits a unique solution in p0, p T q. Moreover, the following properties hold where T will be fixed later on as a function of R. We aim to construct a mapping S : X R,T Ñ X R,T which results to be a contraction for a suitable choice of T , with respect to the norm of Cpr0, T s; L 2 pΩqq X L 2 p0, T ; H 2 pΩqq. As a consequence, we will deduce that it admits a unique fixed point and we will show that this fixed point provides a solution to Problem pP q. By abusing of notation, in the following we will use the same symbol C for possibly different positive constants, depending only on the data of the problem, on Ω and (continuously) on T at most; while CpRq will denote a positive constant which also depends on R. Consequently, by (3.3), we can define an operator S 1 pp sq " p c. Then, exploiting the bounds on c 0 (see (2.5)) and on A and B (see (2.9)), we deduce So that, we have pA`Bc 0 pxqqe A´t 0 p spx,τ qdτ´B c 0 pxq ą A`Bc 0 pxq´Bc 0 pxq ą A ą 0, in Ωˆp0, T q, and then (cf. also (3.3)) As a consequence we deduce, for all px, tq P Ωˆp0, T q, the following inequalities Thus, we eventually obtain 3.2. Second step: finding r. We indicate by p s |Γ the trace on Γ of the fixed p s P X R . Analogously, p c |Γ indicates the trace of p c (given by (3.3)) on Γ. By construction of X T and (3.11), we have that p s |Γ , p c |Γ P L 2 p0, T ; H 3{2 pΓqq. Here, for the sake of simplicity, we use the same notation p c and p s on Γ also for their traces. Owing to standard theory for ordinary differential equations, exploiting (2.2) it is straightforward to find a unique p r " S 2 pp s, p cq solving the Cauchy problem, for almost any x P Γ, Testing (3.12) (replacing r with p r) by p r t and integrating over p0, tq, by means of Young's inequality, the smoothness conditions on Ψ (see (2.2)) and the chain rule for the sub-differential BW , we obtain (see (2.7)) Thus, using the Gronwall lemma and recalling (3.9), we eventually get We first deal with the first integral on the left hand side, integrating by parts in time and exploiting (2.1) and (3.2). So that we obtain t 0ˆΩ pφpp cqsq t s "ˆΩ φpp cptqqs 2 ptq´ˆΩ φpc 0 qs 2 0´ˆt 0ˆΩ φpp cqss t (3.16) Combining (3.15) with (3.16), we get (see (3.8) and (3.11)) from which we deduce and then We can now use a generalized Gronwall's lemma (see [2, Teorema 2.1]). Note that }p c t } L 8 pΩq is bounded in L 2 p0, T q by (3.11). We obtain }s} L 8 p0,T ;L 2 pΩqq ď CpRq, (3.20) and, thanks to (3.18), we get Recalling that s is a strong solution, we test equation (2.15) by s t , where p c " S 1 pp sq and p r " S 2 pp s, S 1 pp sqq. Integrating over p0, tq we obtain Then, on account of (3.23) and (3.24), from (3.22) we deducê Observe now thaťˇˇˇˆt Recalling (3.14) and (3.21) we geťˇˇˇˆt Then, combining (3.26)-(3.30) with (3.25), we find m 4 }s t } 2 which gives (see also (3.20)) Actually, we can now estimate s t in L 8 p0, T ; L 2 pΩqq. To this aim let us (formally) differentiate the equation for s with respect to time and get φpp cqs tt`2 Bp c t s t`B p c tt s´div pφpp cq∇sq t (3.33) "´Bp c t p cs´φpp cqp c t s´φpp cqp cs t .
Observe that we can determine s 1 " s t p0q from the identity (cf. (1.1)) and, due to the regularity of s 0 and c 0 , we can deduce that s 1 P L 2 pΩq. The boundary condition turns out to be rewritten as (see Remark 2.2) In order to make our argument rigorous, we can observe, for instance, that the resulting problem (3.33)-(3.35) (written with respect to the unknowns " s t ) has the same structure of our original problem, and thus it admits a unique solution in the standard energy space. We thus test (3.33) by s t and integrate over p0, tq, obtaining Recalling the properties of p c, p r and s, from (3.37) we obtain the estimate pφpp cqsp cq t ss t . On account of the Young inequality and the Sobolev injection H 1 pΩq ãÑ L 4 pΩq, observe that it holds Combining (3.38) and (3.39) we then find Here, we have used the Young inequality, the trace theorem and the Sobolev injection once more. In addition note that, by (3.11), (3.14), and (3.32), from (3.40) we deduce Then we can apply the generalized Gronwall lemma, getting Assuming p T P p0, T s sufficiently small such that so that S : X R, p T Ñ X R, p T . It is clear from the construction that a fixed point of S is a local (in time) solution to our problem.
3.4. The contraction argument. Let us fix p s i P X R, p T , i " 1, 2, and set p c i " S 1 pp s i q, p r i " S 2 pp s i , p c i q and s i " S 3 pp c i , p r i q " Spp s i q. We first estimate pp c 1´p c 2 q and pp r 1´p r 2 q in suitable norms. To this aim, let us first write (3.2) for i " 1, 2, take the difference and test by p c 1´p c 2 . After integrating over p0, tq and using the Young inequality, we get (see (3. Using once more the generalized Gronwall lemma and the fact that }p s 2 } L 8 pΩq is bounded in L 2 p0, p T q (see (3.46)), we obtain }pp c 1´p c 2 qptq} L 2 pΩq ď CpRq}p s 1´p s 2 } L 2 p0,t;L 2 pΩqq . In order to estimate pp r 1´p r 2 q, we write (3.12) for the two indices, take the difference and test by pp r 1´p r 2 q. Then, we integrate in time and exploit the monotonicity of W 1 and (2.2). This implies Thus, due to (3.58), we can eventually deduce }pp r 1´p r 2 qptq} L 2 pΓq ď CpRq}p s 1´p s 2 } L 2 p0,t;H 1 pΩqq . (3.60) Consider now the two equations for s 1 and s 2 . Subtracting the second from the first, we get pφpp c 1 qs 1 q t´p φpp c 2 qs 2 q t´d iv pφpp c 1 q∇s 1´φ pp c 2 q∇s 2 q (3.61) Recalling the expression for φ, we obtain with the boundary condition φpp c 1 qB n s 1´φ pp c 2 qB n s 2 "´νpp r 1 qps 1´s q`νpp r 2 qps 2´s q. (3.63) Testing (3.62) by s 1´s2 and integrating over p0, tq, we find pp c 1´p c 2 q∇s 2 ∇ps 1´s2 qˆt 0ˆΓ νpp r 1 qps 1´s2 q 2`p νpp r 1 q´νpp r 2 qqps 2´s qps 1´s2 q "´ˆt 0ˆΩ φpp c 1 qp c 1 ps 1´s2 q 2ˆt 0ˆΩ`A pp c 1´p c 2 q`Bpp c 1´p c 2 qpp c 1`p c 2 q˘s 2 ps 1´s2 q.
Observe that

66)
and }pp c 1 q t } L 8 pΩq is bounded in L 2 p0, p T q (see (3.11)). Then, it holds (see (3.11) and (3.32))ˇˇˇBˆt 0ˆΩ pp c 1´p c 2 qps 2 q t ps 1´s2 qˇˇˇˇď Cˆt for a sufficiently small δ to be chosen later. Analogously, we have (see (3.56 The gradient terms are estimated using (3.58), namely,ˇˇˇBˆt 0ˆΩ pp c 1´p c 2 q∇s 2 ∇ps 1´s2 qˇˇˇˇ(3.69) Observe that we can control the boundary term in this wayˇˇˇˆt Combining (3.64)-(3.70) and taking δ sufficiently small, we deduce }ps 1´s2 qptq} 2 L 2 pΩq`ˆt 0 }ps 1´s2 q} 2 Finally, making use of the generalized Gronwall lemma, we get }ps 1´s2 qptq} 2 Thus, for a suitable large j P N (we recall that R is fixed) we have that S j is a contraction in X R, p T . Hence, S admits a unique fixed point s, which is the local in time solution to our original problem. 4. Global existence result for the initial problem. In this section, we show that (2.10) allows us to obtain a global a priori L 8 bound on s. Therefore, the previous fixed point argument can now be carried out without restriction on T so that the unique fixed point is solution to our problem on the whole time interval p0, T q.
Then we can apply once more [14, Lemma 2.1] and deduce that z ě 0 in Ωˆp0, T q. So that (see (3.50)) 0 ď spx, tq ď S 0 , for a.e. px, tq P Ωˆp0, T q. Here C also depends on S 0 (see (2.10)). Then, we can test (2.15) by s and integrate over p0, tq. Arguing as in (3.19)-(3.20) we now deduce, in particular, that where now C does not depend on R.
By means of the previous estimates and the generalized Gronwall lemma we also get that }c} L 8 p0,T ;H 1 pΩqq ď C, (4.9) where C does not depend on R.
We can now argue as for (3.13), where the norms }s} L 2 p0,T ;L 2 pΓqq and }c} L 2 p0,T ;L 2 pΓqq are controlled by the norms in }s} L 2 p0,T ;H 1 pΩqq , }c} L 8 p0,T ;H 1 pΩqq , respectively. Thus, we obtain (replacing p r with r) }r} H 1 p0,T ;L 2 pΓqq ď C, for some C independent of R. Choosing now, in the fixed point argument, the constant R ą 0 large enough such that pT`T 1{2 qC`}s 0 } L 2 pΩq ď R, then we can prove that S : X R,T Ñ X R,T is a contraction in X R,T so that the solution is defined on p0, T q.

Numerical examples.
A fully implicit finite elements scheme has been used to solve numerically the proposed model where s, c and r are defined within the body or along the boundary Γ of the domain Ω. The evolution equation for s (see (1.1)) is discretized with linear finite elements. The code is based on the freely available Open Source package deal.II [3].
Computational tests were performed to verify the efficiency of the numerical technique in 2D setup (x " x 1 e 1`x2 e 2 ). Moreover, these computations permit to illustrate the main features of the proposed approach. The considered domain is a simple marble square domain 1ˆ1 mm 2 invested by a polluted air flow along the left vertical side whereas the other faces are isolated. The distribution of the pollutant SO 2 on the external border is constant with a concentration equal to s 0 . Moreover, the material is homogeneous so that the initial condition for the calcite density reads c px, 0q " c 0 , being c 0 a positive constant, whereas the concentration of SO 2 within the solid is null so s px, 0q " 0.
The numerical examples have been obtained by assuming simple expressions for the functions of (1.3), (1.4). In particular, according to the physical evidence that the higher is the rugosity value more exchange surface is available, two monotone relations have been considered for νprq with r ě 0: ‚ linear ν prq " ν 0`ν l´ν0 r l r ‚ parabolic ν prq " ν 0`ν l´ν0 r l 2 r 2 where the constants are ‚ ν 0 : the minimum value for a fully flat surface (r " 0) ‚ ν l : the maximum value when it is assumed that r P r0, r l s (cf. Remark 2.3, setting r l " R 0 ).
The two curves give the same values of ν prq for r P t0, r l u.
For simplicity, no constraints on the variable r are introduced so to have BW prq " Ψ 1 prq " 0 and environmental sources are not considered thus F " 0. Under these hypotheses, equation (1.4) reduces to r t`G pr, c, sq " 0, on Γˆp0, T q. (5.1) For the function Gpr, c, sq the following expression is adopted Gpr, c, sq "´ϕ pcq csˆ1`r 1`r˙g , being g a parameter that defines the rate of rugosity evolution and has to be fitted according to experimental evidences. In the numerical examples different initial rugosity distributions have been considered: piecewise constant or random distribution along the boundary. In both cases linear and parabolic diffusion functions for νprq have been adopted. Finally, it has been assumed λ " 100, g " 30 and the time is discretized into n constant time intervals ∆t " 1{5000. 5.1. Piecewise constant rugosity. In the first example, piecewise constant value of the rugosity is assumed along the contour: for x 2 ă 0.5, r px, 0q " 0.5r 0 and for x 2 ě 0.5, r px, 0q " 2r 0 with r 0 a positive constant. In Figs. 1, 2, 3 are reported the profiles of the three variables c, s, r along the vertical side invested by the pollutant SO 2 for several time steps. Both linear and parabolic evolution laws for νprq are considered.
The sulphation process reported in Fig. 1, due to the different assumed values of r, evolves with various velocities along the boundary. In fact, for small time step values, the profile of c presents a significant gradient near x 2 " 0.5 where the rugosity is initially discontinuous. This phenomenon is much more evident in case the parabolic relation for νprq is assumed. Subsequently, the solution becomes smoother and tends to be homogenous as c Ñ 0. This behavior can be understood analyzing the evolution of s and r as plotted in Figs. 2, 3 respectively. In fact, small rugosity values act as a barrier to the penetration of s. Even for these two variables high gradients are evident near the middle of the side. Indeed, the parabolic function for νprq reveals higher gradient values. The effect of different rugosity values also influences the solution within the solids. The evolutions of the variables c and s are plotted along two horizontal lines located at x 2 " t0.25, 0.75u for different time steps in Figs. 4, 5. At x 2 " 0.25, where the initial rugosity is high, the sulphation process is affected by the rapid diffusion of s within the solid that activates the calcite transformation. On the contrary, the diffusion of SO 2 at x 2 " 0.75, and consequently the marble degradation, begins slowly due to the initial small value of νprq. For large time steps (n ą 100) the influence of the surface rugosity turns to be negligible on the inner solutions. The maps of s within the solid is plotted in Fig. 6 for time steps n " t5, 15u and parabolic νprq. In the lower portion of the domain, due to higher rugosity values along the border, larger diffusion of SO 2 is evident. 5.2. Random rugosity. Subsequently, a random surface rugosity based on Weibull's statistics is assigned on Γ. The initial value ofr 0 is therefore assumed as with r 0 and m being, respectively, the Weibull shape and modulus parameters, and λ a random variable ranging from 0 to 1. Here, r 0 is the mean value of surface rugosity. In the computation m " 10 has been considered. The initial rugosity profile is reported in Fig. 9 (the profile at timestep n " 0). In the evolution, the asperities are maintained for νprq linear whereas are emphasized in case of parabolic relation. For large time the rugosity does not evolve anymore and the heterogenous profile is preserved.
The initial homogenous distribution of c is abandoned during the sulphation process as clearly stated by Fig. 7 that reports  νprq. Subsequently, the asperities vanish as c Ñ 0. Analogous behavior appears for the SO 2 concentration as reported in Fig. 8. During the sulphation process, the rugosity profile, as stated in Fig. 9, maintains the initial shape with an increment of the fluctuation. 6. Conclusions. The local in time well-posedness and the global existence theorem are preliminary but essential results related to our problem. A nontrivial issue might be the analysis of the longtime behavior of solutions and, in particular, their possible convergence to a steady state. Moreover, it would be interesting (and challenging) to extend our theoretical analysis to significant generalizations which account for further features of this complex phenomenon. In particular, in the coupling of bulk-surface equations, further actions should be taken into account, like the bulk damage and thermo-mechanical effects forcing the surface deterioration. The presented numerical simulations highlight the main features of the proposed model. The introduction of the superficial rugosity permits to account for an important physical property of the material that strongly influences initiation and the evolution of the sulphation stone subjected to pollutant. In fact, high rugosity values offers wider surface for the SO 2 diffusion that induces the sulphation process. On the other hand, extremely smooth surfaces slow down the degradation process.
At the present stage, two natural developments are possible. Firstly, an extensive experimental campaign has to be performed in order to fit the parameters of the proposed model and to define a proper function for νprq. Secondly, the fracture phenomenon and mechanical degradation has to be coupled with the process of sulphation of the stones. In particular, the best candidate is the variational formulation of fracture mechanics that recently has gained much attention and demonstrated to be effective to determine complex crack paths [6].