SIAM News

	Figure 1. The Caryatid Porch of the Erechtheion. From Wikipedia, en.wikipedia.org/wiki/Caryatid. Click to enlarge.

by Antonio Fasano and Roberto Natalini

   The picture of the caryatids in the celebrated Erechtheion on the Acropolis (Figure 1) stands as heartbreaking evidence of an invaluable cultural heritage that we have lost forever.

   By this we mean that the current generation is responsible for the deterioration of the caryatids, as well as for the devastation of many other stone masterpieces that, until recently, had remained intact for thousands of years. Paradoxically, even as we participate in their destruction, we respect these masterpieces from the past and take pride in scientific methods devised to preserve them.

   The list of great monuments and frescoes destroyed or seriously damaged by negligence, or even purposely, has been painfully long for many years. The most lethal attacks, however, have come in the last few decades, from a silent enemy: air pollution. Temples like the Taj Mahal in Agra have come to us with their finest details still intact because they are situated in non-polluted areas.

   Chemistry, of course, has the most to reveal about the causes of marble degradation. A survey of the chemical (and biological) processes involved in both the degradation and the restoration of marble can be found in [4]. But what can mathematics say about the cruel struggle of marble monuments against air pollution?

   Mathematical modeling can lead us to a better understanding of the relative importance of the simultaneous processes involved. Models can help to predict the evolution of the phenomenon, and its dependence on seasonal variations in the basic physical quantities (air pressure, temperature, moisture content, rainfall, and—naturally—the concentrations of the pollutants). As a consequence, a mathematical model can serve as a basis for determining an optimal strategy for restoration or even prevention. Continuing one step further, a model can be used to predict the outcome of specific restoration techniques.

The Silent Killer: SO₂

   The substance most aggressive in attacking calcium carbonate (CaCO₃), the basic ingredient of marble, is sulfur dioxide (SO₂). The reaction is not direct, but rather is mediated by water vapor present in air. The chain of reactions begins with the creation of H2SO₄, which then reacts with CaCO₃, producing calcium sulfate (CaSO₄) and CO₂. A reduced scheme for the reactions (not showing the intermediate steps) is

CaCO₃ + SO₂ + O₂ + 2H₂O                          
→ CaSO₄ · 2H₂O + CO₂ .                         (1)

The precious carbon atom, stolen from the stone, passes into the air—it is in this way that the classical beauty of the caryatids has been dispersed into the atmosphere. An-other similar reaction can take place, yielding calcium sulfite (CaSO₃), but we disregard it here for simplicity.

   The layer of calcium sulfate left in place is called gypsum. The polished surface of the undamaged marble has given way to the opacity of gypsum, but the decline in aesthetic appeal is less important than the mechanical damage incurred: Gypsum is much more porous than marble and, unlike marble, is crumbly and easily eroded by the action of wind and rain. At critical points, moreover, the expansion resulting from the increased porosity can produce stresses large enough to fracture the gypsum.

Figure 2. SEM–EDS map of the chemical elements on a sample surface exposed for 144 hours (right: calcium; left: sulfate), as obtained in a laboratory test in [2]. Click to enlarge.

A Mathematical Approach to Sulfation

   The question of space and time scales is a key consideration in setting up a mathematical model. Because SO₂ is (fortunately) very diluted in the atmosphere—about 0.005 ppm in urban areas, according to the U.S. Environmental Protection Agency—the destruction occurs very slowly, even though the sulfation reaction can be considered instantaneous and proceeds along a sharp front in stone of very low porosity—like very high-quality marble, as observed experimentally in [2,4] (see Figure 2).

   The typical penetration speed of the front is of the order of several microns per year (if the gypsum layer has not been removed). For surfaces in which the curvature is not too great, a one-dimensional model is therefore appropriate.

Transport of SO₂ and Front Advance. Analytical and numerical results established for a preliminary mathematical model are presented in [1,3]. Here, we refine that analysis, using the frame of reference sketched in Figure 3.

   We denote by s the concentration of SO₂ in the pores of the gypsum. It is quite reasonable to assume that the SO₂ flux relative to air in the gypsum pores is Fickian, i.e., that it is driven by the gradient of s. Such a gradient is evidently created by the reaction at the sulfation front, where the consumption rate of SO₂ is proportional to the velocity of the front.

   The flux intensity of SO₂ can then be expressed by

j_S = \varphi_{\gamma} (-d_S \partial_X s + \upsilon_a s), \qquad (2)

where \varphi_{\gamma} is the porosity of gypsum (the volume fraction occupied by the pores), which for simplicity we consider constant; d_S is the diffusivity of SO₂ in air, and \upsilon_a is the average velocity of air molecules. The equation governing SO₂ transport is thus

\varphi_{\gamma} \partial_t s + \partial_X j_s = 0. \qquad (3)

   On the reaction front x = \sigma(t), SO₂ reacts completely, because the reaction is instantaneous, provided that enough water is available. Indeed, here we consider an infinite reaction rate, which, according to results in [3], is an excellent approximation. We can thus say that

s (\sigma(t), t) = 0. \qquad (4)

   Balancing mass on the reaction front yields

-\varphi_{\gamma} \frac{d_s}{M-S} \partial_x s = \frac{\rho_m}{M_m} \dot{\sigma} , \qquad (5)

where M_s and M_m are the molar weights of SO₂ and CaCO₃, respectively, and \rho_m is the density of pristine marble. In writing (5), we simply counted the molecules involved in the reaction: On the left we have the number of moles of SO₂ reaching the unit surface of the front per unit time; on the right is the corresponding number of moles of CaCO₃ that disappear. (The molar ratio of the two species in the reaction is 1:1, as shown by (1).) We return later to the question of how (4) and (5) must be modified when the process is slowed by moisture shortages.

   The pore concentration of SO₂ at the outer surface of the gypsum layer x = \sigma_0(t) equals the concentration of SO₂ in air:

s(\sigma_0 (t), t) = s_a \qquad (6)

Motion of Gypsum. Because the ratio \omega between the molar densities of marble (\rho_m /M_m) and gypsum (\rho_{\gamma}/ \rho_{\gamma}) greater than 1, the gypsum layer possesses a velocity:

\dot{\sigma}_0 = -\omega \dot{\sigma}. \qquad (7)

 

	Figure 3. Growth of the internal and external free boundaries. Click to enlarge.

Thus, if \sigma_0(0) = \sigma(0) = 0, the outer surface obeys

\sigma_0(t) = -\omega \sigma(t), \qquad (8)

and the thickness of the gypsum layer is

\sigma(t) - \sigma_0(t) = (1 + \omega) \sigma(t), \qquad (9)

   Because \omega – 1 is not small, this swelling effect cannot be neglected.

Motion of Air. In principle, the mechanism responsible for the motion of air within the gypsum is rather complicated, because the sudden expansion of the medium following the reaction creates a pressure drop at the front (the amount of air present in the pristine marble being negligible). A pressure gradient thus arises, forcing air through the gypsum to fill the volume made available. The flow of air relative to gypsum can be described by Darcy’s law. Rescaling with the appropriate reference quantities, however, makes it clear that in most cases, at the space and time scales of interest, we can assume the air density throughout the growing gypsum layer to be basically uniform and constant (we neglect changes in external conditions, like the action of wind). In short, there is no relative motion between air and gypsum at any time, so that

\upsilon_a \simeq \dot{\sigma}_0. \qquad (10)

Role and Motion of Water Vapor. In considering the role of water, we begin with a question: Why must the relative humidity of air exceed a given threshold (usually fixed by practitioners at 75%) for the sulfation process to occur?

   As shown in (1), conversion of a molecule of CaCO₃ into a molecule of CaSO₄ requires that one molecule of SO₂ and two molecules of H₂O simultaneously reach the molecule of CaCO₃. Because H₂O is diluted, and SO₂ extremely diluted, in air, the probability of such a triple encounter at the reaction front is very small. The implication is that for sulfation to occur, water must be constantl  present on the marble surface as a liquid film, so that the arrival of any molecule of SO₂ at the front will trigger the reaction.

   Indeed, marble is hygroscopic, and the creation of a water film takes place through a sorption–desorption mechanism. Depen-ding on the balance between sorption and desorption (which is clearly influenced by temperature and relative humidity), the film may or may not have the continuity and thickness required to keep the reaction going as long as SO₂ is available, even at small concentrations. The rate constants of sorption–desorption are much larger than the reciprocal time scale of the overall sulfation process; the process of building up the water film can thus be considered independent of sulfation, as long as the air is not too dry. We consider this a reasonable explanation for the moisture threshold needed for marble sulfation. When the air is dry enough to preclude any moisture coating, of course, no sulfation can occur and \dot{\sigma} = 0, implying that \partial s / \partial x = 0 on x = \omega, while s = 0 no longer applies. When the relative humidity is close to the threshold, an intermediate behavior is observed  The water film breaks down into humid spots (which have their own dynamics), producing on average a slowing of the front. The corresponding free boundary conditions, replacing (4) and (5), can be written as follows:

\frac{j_s}{M_s} = \frac{\rho_m}{M_m} \dot{\sigma} + \varphi_{\gamma} \frac{s}{M_s} \dot{\sigma},

\frac{\rho_m}{M_m} \dot{\sigma} = \varphi_{\gamma} \alpha \frac{s}{M_s}, \qquad (11)

where \alpha measures the efficiency of the reaction and depends on the concentration of H₂O in air. The first equation expresses the total molar balance of SO₂ (the total flux accounts for the mass loss rate in the reaction and for the advection of the residual SO₂ by the moving front). The second condition specifies the reaction speed. The coefficient \alpha tends to +\infty when the relative humidity approaches the full-speed threshold, and to 0 when the relative humidity decreases to the lower threshold of no reaction. In the first case the concentration s on the free boundary is forced to go to 0 and we recover the pair (4), (5) characterizing the full-speed regime. With \alpha = 0, by contrast, the front stops and at the same time the SO₂ flux vanishes, yielding \partial x s = 0.

   For the transport of H₂O, using w to denote the pore concentration, we can write the total flux as

j_w = \varphi_{\gamma}(-d_w \partial_x w + w \upsilon_a), \qquad (12)

with d_w denoting the diffusivity of H₂O in air, and the balance equation as

\varphi_{\gamma} \partial_t w + \partial_x j_w = 0. \qquad (13)

   When sulfation proceeds at full speed and the water film is steady (i.e., the relative humidity is above the threshold and temperature is constant), the counting of H₂O molecules at the reaction front gives

\frac{j_w}{M_w} = 2 \frac{\rho_m}{M_m} \dot{\sigma} + \varphi_{\gamma} \frac{w}{M_w} \dot{\sigma},  \qquad (14)

where M_w is the molar weight of water; the factor 2 indicates that two molecules of H₂O enter the reaction.

   Equation (14) has the same meaning as the first equation of (11), and it is valid for any regime. In (14), because of the very low porosity of pristine marble, we neglect any water molecules that might be present in the marble. The value of w at the outer surface x = \sigma_0(t) equals the concentration w_a of H₂O in air.

   Clearly, the introduction of temperature variations, with the consequent dynamics in the liquid water film coating the reaction front, would make the H₂O balance at the interface considerably more complicated.

Conclusions and Perspectives

   As discussed in this article, the sulfation of marble involves not only the chemical reaction at the front, but also nontrivial transport phenomena for SO₂, H₂O, and air (which we neglected completely). We also neglected the destiny of CO₂ (one of the reaction products), which may enter into other reactions, and, above all, the parallel actions of other impurities, like NO₂, present with SO₂ in polluted air. Again, we refer readers interested in the complete chemical picture to [4].

   The story also has a biochemical side, which is particularly interesting for mathematicians. Indeed (again, see [4]), one of the techniques used in marble remediation is the implantation of a population of bacteria (Desulfovibrio desulfuricans) capable of metabolizing CaSO₄; in effect, the bacteria set in motion the inverse reaction by which CaSO₄ reverts to CaCO₃: a perfect subject for experts in population dynamics.

   Despite the complexity of the model sketched here, there may be cases in which, after rescaling, the problem can be greatly simplified. If, for instance, (10) is valid and H₂O is present in large enough quantities, and the physical constants are such that we can neglect the inertia term in (3), then j_s is a function only of time and its value can be deduced from (5):

\varphi_{\gamma} (d_s \partial_x s + \omega \dot{\sigma} s )= \frac{M_s}{M_m} \rho_m \dot{\sigma}. \qquad (15)

   Integrating (15) with respect to x with the condition (6) (e.g., with s_₀ constant), we arrive at the expression

s = -\sum + (s_0 + \sum)e^{-\omega \dot{\sigma}\frac{x - \sigma_0}{d_s}},

with \sum = \frac{M_s \rho_m}{M_m \varphi_{\gamma} \omega}, \qquad (16)

Imposing s(\sigma(t), t) = 0 and recalling (9), we deduce a simple ordinary differential equation for \sigma:

\omega (1 + \omega) \frac{\sigma \dot{\sigma}}{d} = \log(1 + \frac{s_0}{\sum}), \qquad (17)

to be integrated with the condition \sigma(0) = 0:

\sigma(t) = \displaystyle \left| \frac{2d_s t}{\omega (1 + \omega)} \log(1 + \frac{s_0}{\sum})\right|^{1/2}

\simeq \displaystyle \left( \frac{2d_s t}{\omega (1 + \omega)} \frac{s_0}{\sum} \right)^{1/2}, \qquad (18)

considering that s_0 / \sum << 1. The formula is slightly different if s_0 depends on t:

\sigma(t) = \left[ \frac{2d_s}{\omega (1 + \omega)} \frac{1}{\sum} \int_0^t s_o(\tau) d\tau \right]^{1/2}, \qquad (19)

   The t^1/2 behavior is in sharp agreement with experimental results in [2,4] and with the numerical simulations in [1], and can be attributed to the Stefan-like nature of the free boundary conditions (4),(5), which is also seen in other phenomena (e.g., the penetration of dry ground by water, i.e., the Green–Ampt problem, well known to soil engineers).

   Equation (18) gives a very expressive answer, emphasizing the aggressiveness of SO2 against a fresh marble surface, but also showing the shielding effect of the gypsum layer. This might suggest the optimal strategy for restoration by gypsum removal: If repeated t^1/2-type kinks in the sulfation front are to be avoided, the gypsum removal must not be done too frequently.

   Another useful suggestion from (18) concerns the choice of space and time scales. If s^*₀ = 14.3 μg/m³ = 0.005 ppm is, say, the average SO₂ concentration on a day in an urban area in the U.S. (averaged EPA data for 2001) and t* = 365 days, a suitable length scale is

\sigma^* = \left| \frac{2d_s t^*}{\omega (1 + \omega)} \frac{s^*_0}{\sum} \right|^{1/2} \simeq 23.6 \mu m;

using the data d_s = 0.1119 ×10^-4 m²/sec,  \varphi_{\gamma} = 30%,

\sum = \frac{M_s \rho_m}{M_m \varphi_{\gamma} \omega} = 3 \mu g / m^3,

\omega = \frac{\rho_m M_{\gamma}}{M_m \rho_{\gamma}} = 2.

With the rescaled variables \bar{s}_0 = s_0 / s^*_0, \bar{t} = t/t^*, \bar{\sigma} = \sigma / \sigma^*, (18) becomes

\bar{\sigma}(\bar{t}) = \sqrt{\tilde{t}\tilde{s}_0}, \qquad (20)

and the same quantities can be used to re-scale the general problem. According to (20), removing the gypsum crust after n years of exposure to the SO₂ concentration s^*₀ causes the loss of n^1/2 × 23.6 μm, while the total loss with yearly restoration over the same period is n × 23.6 μm. Thus, the damage increases by a factor of 5 in only 25 years! But the length scale σ^* depends on (s^*₀)^1/2; to halve the length σ^* of the crust after one year, therefore, we need to decrease the SO₂ concentration to a quarter of the initial averaged value.

   We consider this a fascinating research area in which mathematicians can produce not only elegant theories, but also very concrete answers. Of course, mathematics cannot console us for the loss of the delicate beauty of the caryatids, but it might help to preserve as yet undamaged treasures for the joy of future generations.

References

[1] D. Aregba-Driollet, F. Diele, and R. Natalini, A mathematical model for the SO2 aggression to calcium carbonate stones: Numerical approximation and asymptotic analysis, SIAM J. Appl. Math., 64 (2004), 1636–1667.

[2] E. Borrelli, C. Giavarini, M. Incitti, M. L. Santarelli, and R. Natalini, A material model for the evolution of gypsum crusts: Nu-merical and experimental results, Proceedings of the 10th International Congress on Deterioration and Conservation of Stone, D. Kwiatkowski and R. Löfvendahl, eds., Stockholm, Sweden, I (2004), 35–42.

[3] F.R. Guarguaglini and R. Natalini, Fast reaction limit and large time behavior of solutions to a nonlinear model of sulphation phenomena, Commun. Partial Differential Equations, to appear.

[4] K. Lal Gauri and J.K. Bandyopadhyay, Carbonate Stone: Chemical Behavior, Durability, and Conservation, John Wiley & Sons, Hoboken, NJ, 1999.

Antonio Fasano is a professor of analytical mechanics at the University of Florence. Roberto Natalini is a mathematician at the Istituto per le Applicazioni del Calcolo of CNR, in Rome.

« Back to July/August 2006 Index