A model for lime consolidation of porous solids^††thanks: This work was supported by the GAČR Grant No. 20-14736S; the European Regional Development Fund Project No. CZ.02.1.01/0.0/0.0/16_019/0000778; and the Austrian Science Fund (FWF) Project V662.

In this paper we investigate mathematically a process which is used by the building industry in order to protect and conserve cultural goods and other structure works. Such structures which are subject to weathering can be strengthened by filling the pores by a water-lime-mixture. The mixture penetrates into the pore structure of the stone and the calcium hydroxide reacts with carbon dioxide and builds calcium carbonate and water. A solid layer is built in the pore space which strengthens the material. The main problem is, however, that the consolidated layer is in practice rather thin and is located too close to the active boundary.

In order to avoid possible ambiguity, it is necessary to explain that the term ‘consolidation’ is to be interpreted here as a process of formation of calcium carbonate which has the property of binding the particles together. It might also be called ‘cementation’ or ‘compaction’.

In the literature there are a lot of works dealing with those problems. Different practical strategies for the wetting and drying regime which lead to a more uniform distribution of the consolidant are compared in [25]. A mathematical model is proposed in [31, 32], and [13], where the authors derive governing equations for moisture, heat, and air flow through concrete. A numerical procedure based on the finite element method is developed there to solve the set of equations and to investigate the influence of relative humidity and temperature. It is shown that the amount of calcium carbonate formed in a unit of time depends on the degree of carbonation, i. e., the availability of calcium hydroxide, the temperature, the carbon dioxide concentration and the relative humidity in the pore structure of the concrete. An extension of the aforementioned papers by studying the hygro-thermal behavior of concrete in the special situation of high temperatures can be found in [17].

In the present case chemical reactions take place. Various approaches exist which describe such processes by models which stem from different backgrounds (e. g., from mixture theory or empirical models). An overview on the development of theories especially for porous media including chemical reactions is given in [14]. The interactions between the constituents of a porous medium are not necessarily of chemical nature which would lead to a chemical transformation of one set of chemical substances to another. Simpler is the mass exchange between the constituents by physical processes like adsorption. Adsorption-diffusion processes have been studied by B. Albers (the former name of B. Detmann) e. g. in [2]. Other works on sorption in porous solids including molecular condensation are [5] or [6]. In these works the diffusivities of water and carbon dioxide are assumed to be strongly dependent on pore humidity, temperature and also on the degree of hydration of concrete. The authors realized that the porosity becomes non-uniform in time. This is an observation which is interesting also in the present case because the structure of the channels clearly changes with the progress of the reaction. A survey of consolidation techniques for historical materials is published in [27]. The influence of the particle size on the efficiency of the consolidation process is investigated in [35]. Experimental determination of the penetration depth is the subject of [7, 8, 9]. Different variants of the consolidants is studied in [15, 22, 28, 33], and an experimental work on mechanical interaction between the consolidant and the matrix material is carried out in [21].

A further work dealing with chemical reactions and diffusion in concrete based on the mixture theory for fluids introduced by Truesdell and coworkers is by A. J. Vromans et al. [37]. The model describes the corrosion of concrete with sulfuric acid which means a transformation of slaked lime and sulfuric acid into gypsum releasing water. It is a similar reaction we are looking at. A similar topic is dealt with in [36], where it is shown how the carbonation process in lime mortar is influenced by the diffusion of carbon dioxide into the mortar pore system by the kinetics of the lime carbonation reaction and by the drying and wetting process in the mortar.

Experimental results of $\mathrm{CaCO_{3}}$ precipitation kinetics can be found in [30]. The porosity changes during the reaction. This was studied by Houst and Wittmann, who also investigated the influence of the water content on the diffusivity of $\mathrm{CO_{2}}$ and $\mathrm{O_{2}}$ through hydrated cement paste [19]. An investigation of the physico-chemical characteristics of ancient mortars with comparison to a reaction-diffusion model by Zouridakis et al. is presented in [38]. A slightly different reaction involving also sulfur is mathematically studied in [10] by Böhm et al. There the corrosion in a sewer pipe is modeled as a moving-boundary system. A strategy for predicting the penetration of carbonation reaction fronts in concrete was proposed by Muntean et al. in [23]. A simple 1D mathematical model for the treatment of sandstone neglecting the effects of chemical reactions is proposed in [12] and further refined in [11].

We model the consolidation process as a convection-diffusion system coupled with chemical reaction in a 3D porous solid. The physical observation that only water can be evacuated from the porous body, while lime remains inside, requires a nonstandard boundary condition on the active part of the boundary. We choose a simple one-sided condition for the lime exchange between the interior and the exterior. The main result of the paper consists in proving rigorously that the resulting initial-boundary value problem for the PDE system in 3D has a solution satisfying natural physical constraints, including the boundedness of the concentrations proved by means of time discrete Moser iterations. We also show the result of numerical simulation in a simplified 1D situation.

The structure of the present paper is the following. In Section 1, we explain the modeling hypotheses and derive the corresponding system of balance equations with nonlinear boundary conditions. In Section 2, we give a rigorous formulation of the initial-boundary value problem, specify the mathematical hypotheses, and state the main result in Theorem 2.2. The solution is constructed by a time-discretization scheme proposed in Section 3. The estimates independent of the time step size derived for this time-discrete system constitute the substantial step in the proof of Theorem 2.2, which is obtained in Section 4 by passing to the limit as the time step tends to zero. A numerical test for a reduced 1D system is carried out in Section 5 to illustrate a qualitative agreement of the mathematical result with experimental observations.

We imagine a porous medium (sandstone, for example) the structure of which is to be strengthened by letting calcium hydroxide particles driven by water flow penetrate into the pores. In contact with the air present in the pores, the calcium hydroxide reacts with the carbon dioxide contained in the air and produces a precipitate (calcium carbonate) which is not water-soluble, remains in the pores, and glues the sandstone particles together. Unlike, e. g., in [18, 34], we do not consider the porosity as one of the state variables. The porosity evolution law is replaced with the assumption that the permeability decreases as a result of the calcium carbonate deposit in the pores. The chemical reaction is assumed to be irreversible and we write it as $Ca(OH)_{2}+CO_{2}\to CaCO_{3}+H_{2}O$.

Notation:

• $\dot{c}^{W}$ … mass source rate of $H_{2}O$ produced by the chemical reaction

• $\dot{c}^{H}$ … mass source rate of $Ca(OH)_{2}$ produced by the chemical reaction

• $\dot{c}^{P}$ … mass source rate of $CaCO_{3}$ produced by the chemical reaction

• $\dot{c}^{G}$ … mass source rate of $CO_{2}$ produced by the chemical reaction

• $m^{W}$ … molar mass of $H_{2}O$

• $m^{H}$ … molar mass of $Ca(OH)_{2}$

• $m^{P}$ … molar mass of $CaCO_{3}$

• $m^{G}$ … molar mass of $CO_{2}$

• $\rho^{W}$ … mass density of $H_{2}O$

• $\rho^{H}$ … mass density of $Ca(OH)_{2}$

• $p$ … capillary pressure

• $s$ … water volume saturation

• ${h}$ … relative concentration of $Ca(OH)_{2}$

• $p^{\partial\Omega}$ … outer pressure

• $h^{\partial\Omega}$ … outer concentration of $Ca(OH)_{2}$

• $v$ … transport velocity vector

• $k(c^{P})$ … permeability of the porous solid

• $q$ … liquid mass flux

• $q^{H}$ … mass flux of $Ca(OH)_{2}$

• $\gamma$ … speed of the chemical reaction

• $\kappa$ … diffusivity of $Ca(OH)_{2}$

• $n$ … unit outward normal vector

• $\sigma(x)$ … transport velocity interaction kernel

• $\alpha(x)$ … boundary permeability for water

• $\beta(x)$ … boundary permeability for the inflow of $Ca(OH)_{2}$

Mass balance of the chemical reaction:

$$\frac{\dot{c}^{P}}{m^{P}}=\frac{\dot{c}^{W}}{m^{W}}=-\frac{\dot{c}^{H}}{m^{H}}=-\frac{\dot{c}^{G}}{m^{G}}.$$

Water mass balance in an arbitrary subdomain $V$ of the porous body:

$$\frac{\,\mathrm{d}}{\,\mathrm{d}t}\int_{V}\rho^{W}s\,\mathrm{d}x+\int_{\partial V}q\cdot n\,\mathrm{d}S=\int_{V}\dot{c}^{W}\,\mathrm{d}x.$$

Calcium hydroxide mass balance in an arbitrary subdomain $V$ of the porous body:

$$\frac{\,\mathrm{d}}{\,\mathrm{d}t}\int_{V}\rho^{H}{h}\,\mathrm{d}x+\int_{\partial V}q^{H}\cdot n\,\mathrm{d}S=\int_{V}\dot{c}^{H}\,\mathrm{d}x.$$

Water mass balance in differential form:

$$\rho^{W}\dot{s}+\mathrm{\,div\,}q=\dot{c}^{W}.$$

Calcium hydroxide mass balance in differential form:

$$\rho^{H}\dot{h}+\mathrm{\,div\,}q^{H}=\dot{c}^{H}.$$

The water mass flux is assumed to obey the Darcy law:

$$q=-k(c^{P})\nabla p,$$

with permeability coefficient $k(c^{P})$ which is assumed to decrease as the amount $c^{P}$ of $CaCO_{3}$ given by the formula

$$c^{P}(x,t)=\int_{0}^{t}\dot{c}^{P}(x,t^{\prime})\,\mathrm{d}t^{\prime}$$

increases and fills the pores.

The flux of $Ca(OH)_{2}$ consists of transport and diffusion terms:

$$q^{H}=\rho^{H}{h}v-\kappa s\nabla{h}.$$

The mobility coefficient $\kappa s$ in the diffusion term is assumed to be proportional to $s$: If there is no water in the pores, no diffusion takes place.

We assume that the transport of $Ca(OH)_{2}$ at the point $x\in\Omega$ is driven by the water flux in a small neighborhood of $x$. In mathematical terms, we assume that there exists a nonnegative function $\sigma$ with support in a small neighborhood of the origin such that the transport velocity $v$ can be defined as

$$v(x,t)=\frac{1}{\rho^{H}}\int_{\Omega}\sigma(x-y)q(y,t)\,\mathrm{d}y.$$

The main reason for this assumption is a mathematical one. The strong nonlinear coupling between $s$ and $h$ makes it difficult to control the bounds for the unknowns in the approximation scheme. We believe that such a regularization of the transport velocity makes physically sense as well.

The wetting-dewetting curve is described by an increasing function $f$:

$$p=f(s).$$

We focus on modeling the chemical reactions. Capillary hysteresis, deformations of the solid matrix, and thermal effects are therefore neglected here. We plan to include them following the ideas of [3] in a subsequent study.

The dynamics of the chemical reaction is modeled according to the so-called law of mass action, which states that the rate of the chemical reaction is directly proportional to the product of the concentrations of the reactants. We assume it in the form

$$\dot{c}^{P}=\gamma m^{P}{h}s(1-s).$$

(1.1)

Its meaning is that no reaction can take place if either no $Ca(OH)_{2}$ is available (that is, $h=0$), or no water is available (that is, $s=0$), or no $CO_{2}$ is available (that is, $s=1$), according to the hypothesis that the chemical reaction takes dominantly place on the contact between water and air. In order to reduce the complexity of the problem, we assume directly that the available quantity of $CO_{2}$ is proportional to the air content.

On the boundary $\partial\Omega$ we prescribe the normal fluxes. For the normal component of $q$, we assume that it is proportional to the difference between the pressures $p$ inside and $p^{\partial\Omega}$ outside the body. For the flux of $Ca(OH)_{2}$, we assume that it can point only inward proportionally to the difference of concentrations and to $s$, and no outward flux is possible. Inward flux takes place only if the outer concentration $h^{\partial\Omega}$ is bigger than the inner concentration $h$:

$$\left.\begin{array}[]{rcl}q\cdot n&=&\alpha(x)(p-p^{\partial\Omega})\\[4.2679pt] q^{H}\cdot n&=&-\beta(x)s(h^{\partial\Omega}-h)^{+}\end{array}\right\}\text{on }\partial\Omega.$$

Let $\Omega\subset\mathbb{R}^{3}$ be a bounded Lipschitzian domain. We consider the Hilbert triplet $V\subset H\equiv H^{\prime}\subset V^{\prime}$ with compact embeddings and with $H=L^{2}(\Omega)$, $V=W^{1,2}(\Omega)$. For two unknown functions $s(x,t),h(x,t)$ defined for $(x,t)\in\Omega\times(0,T)$ the resulting PDE system reads

		$\displaystyle\int_{\Omega}(\rho^{W}s_{t}\phi(x)+k(c^{P})f^{\prime}(s)\nabla s\cdot\nabla\phi(x))\,\mathrm{d}x+\int_{\partial\Omega}\alpha(x)(f(s)-f(s^{\partial\Omega}))\phi(x)\,\mathrm{d}S(x)$
		$\displaystyle\qquad=\gamma m^{W}\int_{\Omega}{h}s(1-s)\phi(x)\,\mathrm{d}x,$		(2.1)
		$\displaystyle\int_{\Omega}(\rho^{H}h_{t}\psi(x)+(\kappa s\nabla{h}-\rho^{H}hv)\cdot\nabla\psi(x))\,\mathrm{d}x-\int_{\partial\Omega}\beta(x)s(h^{\partial\Omega}-h)^{+}\psi(x)\,\mathrm{d}S(x)$
		$\displaystyle\qquad=-\gamma m^{H}\int_{\Omega}{h}s(1-s)\psi(x)\,\mathrm{d}x.$		(2.2)

for all test functions $\phi,\psi\in V$, where $s^{\partial\Omega}:=f^{-1}(p^{\partial\Omega})$, and with initial conditions

$$s(x,0)=s^{0}(x),\quad h(x,0)=h^{0}(x)\quad\mathrm{for\ }x\in\Omega.$$

(2.3)

The meaning of Hypothesis 2.1 (vi) is that the boundary $\partial\Omega$ is inhomogeneous, with different permeabilities at different parts of the boundary. The transport of water (supply of $Ca(OH)_{2}$) through the boundary takes place only on parts of $\partial\Omega$ where $\alpha>0$ ($\beta>0$, respectively).

The remaining sections are devoted to the proof of the following result.

We omit the positive physical constants which are not relevant for the analysis. The strategy of the proof is based on choosing a cut-off parameter $R>0$, replacing $h$ in the nonlinear terms with $Q_{R}(h)=\min\{h^{+},R\}$, $1-s$ with $Q_{R}(1-s)$, and $v$ in (2.2) with $v^{R}:=(Q_{R}(|v|)/|v|)v$. We also extend the values of the function $f$ outside the interval $[0,1]$ by introducing the function $\tilde{f}$ by the formula

$$\tilde{f}(s)=\left\{\begin{array}[]{ll}f(0)+f^{\prime}(0)s&\mathrm{for\ }s<0,\\ f(s)&\mathrm{for\ }s\in[0,1],\\ f(1)+f^{\prime}(1)(s-1)&\mathrm{for\ }s>1,\end{array}\right.$$

and consider the system

		$\displaystyle\int_{\Omega}(s_{t}\phi(x)+k(c^{P})\tilde{f}^{\prime}(s)\nabla s\cdot\nabla\phi(x))\,\mathrm{d}x+\int_{\partial\Omega}\alpha(x)(\tilde{f}(s)-\tilde{f}(s^{\partial\Omega}))\phi(x)\,\mathrm{d}S(x)$
		$\displaystyle\qquad=\int_{\Omega}Q_{R}(h)sQ_{R}(1-s)\phi(x)\,\mathrm{d}x,$		(2.4)
		$\displaystyle\int_{\Omega}(h_{t}\psi(x)+(s\nabla{h}-hv^{R})\cdot\nabla\psi(x))\,\mathrm{d}x-\int_{\partial\Omega}\beta(x)s(h^{\partial\Omega}-h)^{+}\psi(x)\,\mathrm{d}S(x)$
		$\displaystyle\qquad=-\int_{\Omega}{h}s(1-s)\psi(x)\,\mathrm{d}x$		(2.5)

for all $\phi,\psi\in V$. We first construct and solve in Section 3 a time-discrete approximating system of (2.4)–(2.5), and derive estimates independent of the time step. In Section 4, we let the time step tend to $0$ and prove that the limit is a solution $(s,h)$ to (2.4)–(2.5). We also prove that this solution has the property that $s\in[0,1]$, $h$ is positive and bounded, and $v$ is bounded, so that for $R$ sufficiently large, the truncations are never active and the solution thus satisfies (2.1)–(2.2) as well.

For proving Theorem 2.2, we first choose $n\in\mathbb{N}$ and replace (2.4)–(2.5) with the following time-discrete system with time step $\tau=\frac{T}{n}$:

		$\displaystyle\int_{\Omega}\left(\frac{1}{\tau}(s_{j}-s_{j-1})\phi(x)+k(c^{P}_{j-1})\tilde{f}^{\prime}(s_{j})\nabla s_{j}\cdot\nabla\phi(x)\right)\,\mathrm{d}x+\int_{\partial\Omega}\alpha(x)(\tilde{f}(s_{j})-\tilde{f}(s^{\partial\Omega}_{j}))\phi(x)\,\mathrm{d}S(x)$
		$\displaystyle\qquad=\int_{\Omega}Q_{R}(h_{j-1})s_{j}Q_{R}(1-s_{j})\phi(x)\,\mathrm{d}x,$		(3.1)
		$\displaystyle\int_{\Omega}\left(\frac{1}{\tau}(h_{j}{-}h_{j-1})\psi(x)+(s_{j-1}\nabla{h_{j}}{-}h_{j}v^{R}_{j-1})\cdot\nabla\psi(x)\right)\,\mathrm{d}x-\int_{\partial\Omega}\beta(x)s_{j-1}(h^{\partial\Omega}_{j}-h_{j})^{+}\psi(x)\,\mathrm{d}S(x)$
		$\displaystyle\qquad=-\int_{\Omega}{h_{j}}s_{j-1}(1-s_{j-1})\psi(x)\,\mathrm{d}x,$		(3.2)

for $j=1,\dots,n$ with initial conditions $s_{0}=s^{0}$, $h_{0}=h^{0}$, and with $c^{P}_{i}=0$ for $i\leq 0$, with

	$\displaystyle q_{j}(x)$	$\displaystyle=-k(c^{P}_{j}(x))\nabla\tilde{f}(s_{j}(x))$	$\displaystyle\mathrm{for\ }\ x\in\Omega,$		(3.3)
	$\displaystyle v_{j}(x)$	$\displaystyle=\int_{\Omega}\sigma(x-y)q_{j}(y)\,\mathrm{d}y$	$\displaystyle\mathrm{for\ }\ x\in\Omega,$		(3.4)
	$\displaystyle s^{\partial\Omega}_{j}(x)$	$\displaystyle=\frac{1}{\tau}\int_{t_{j-1}}^{t_{j}}s^{\partial\Omega}(x,t)\,\mathrm{d}t$	$\displaystyle\mathrm{for\ }\ x\in\partial\Omega,$		(3.5)
	$\displaystyle h^{\partial\Omega}_{j}(x)$	$\displaystyle=\frac{1}{\tau}\int_{t_{j-1}}^{t_{j}}h^{\partial\Omega}(x,t)\,\mathrm{d}t$	$\displaystyle\mathrm{for\ }\ x\in\partial\Omega,$		(3.6)

where $t_{j}=j\tau$ for $j=0,1,\dots,n$. Moreover, we define inductively

$$c^{P}_{j}-c^{P}_{j-1}=\tau h_{j}s_{j}(1-s_{j})\ \mbox{ for }\ j=1,\dots,n.$$

(3.7)

We now prove the existence of solutions to (3.1)–(3.2) and derive a series of estimates which will allow us to pass to the limit as $n\to\infty$. We denote by $C$ any positive constant depending possibly on the data and independent of $n$.

For $\varepsilon>0$ we denote by $H_{\varepsilon}:\mathbb{R}\to\mathbb{R}$ the function

$$H_{\varepsilon}(r)=\left\{\begin{array}[]{ll}0&\mathrm{for\ }r\leq 0,\\[5.69054pt] \frac{r}{\varepsilon}&\mathrm{for\ }r\in(0,\varepsilon),\\[5.69054pt] 1&\mathrm{for\ }r\geq\varepsilon\end{array}\right.$$

(3.8)

as a Lipschitz continuous regularization of the Heaviside function, and by $\hat{H}_{\varepsilon}$ its antiderivative

$$\hat{H}_{\varepsilon}(r)=\left\{\begin{array}[]{ll}0&\mathrm{for\ }r\leq 0,\\[5.69054pt] \frac{r^{2}}{2\varepsilon}&\mathrm{for\ }r\in(0,\varepsilon),\\[5.69054pt] r-\frac{\varepsilon}{2}&\mathrm{for\ }r\geq\varepsilon\end{array}\right.$$

(3.9)

as a continuously differentiable approximation of the “positive part” function. Note that we have $rH_{\varepsilon}(r)\leq 2\hat{H}_{\varepsilon}(r)$ for all $r\in\mathbb{R}$.

The main issue will be a uniform upper bound for $h_{j}$ which will be obtained by a time-discrete variant of the Moser-Alikakos iteration technique presented in [4]. We start from some preliminary integral estimates of $s_{j}$.