1 Introduction

The partial diffusion differential equation for the uniaxial case is

$$\frac{\partial C}{\partial t} = \frac{\partial }{\partial z}\left[ {D(C)\frac{\partial C}{\partial z}} \right]$$
(1)

Here C is the water concentration, t the time, z the depth coordinate, and D the diffusivity that may depend on the water concentration.

For water vapour as the environment, the surface condition is

$$\frac{dC}{dz} = \frac{h}{D}(C - C_{0} )\;\;{\text{at}}\; z = 0,$$
(2)

where C0 is the concentration of molecular water reached at z = 0 for t → ∞.

Following the suggestion by Doremus ([1], Section 4.7), the phenomenological parameter h in (2) may be interpreted as a reaction parameter for a slow surface reaction that limits the entrance of molecular water species.

On the other hand, a simpler phenomenological description is possible by assuming that a barrier exists to the transport of water across the surface of the glass. The barrier gives rise to a mass transfer coefficient for diffusion, which slows the passage of water into the glass [1]. Each of the assumptions yields the same set of mathematical equations.

The Eqs. (1) and (2) can be solved numerically as was done in [2]. Disadvantage of numerical computations is the fact that the used extend of the z-range must be finite and, consequently, the semi-infinite body can only be approximated. In the following considerations we will give exact and approximate solutions of the diffusion differential Eq. (1) for the half-space.

The problems mentioned here can be solved by “classical” approaches. For much more complicated problems, more specific methods are needed. For example, diffusion is treated in disordered materials by Giona and Roman in [3]. In this reference a solution is analytically obtained for diffusion on fractals proposed for materials showing anomalous transport behavior.

The reason for the concentration-dependent diffusivity used here is the hydroxyl generation at the silica surface exposed to humid environments. Water penetrated into silica reacts with the silica network according to [4]

$${{\equiv}{\text{Si}}} {-} {\text{O}} {-} {{\text{Si}}{\equiv}} + {\text{H}}_{ 2} {\text{O}} \leftrightarrow {{\equiv}{\text{SiOH}}} + {{\text{HOSi}}{\equiv}}$$

with the concentration of the hydroxyl S = [≡SiOH] and that of the molecular water C = [H2O].

A swelling effect in water-containing silica at high temperatures was early reported by Brückner [5, 6], Shackelford [7] and Shelby [8]. These authors showed that the density decreases by the silica/water-reaction. This decrease is equivalent to a volume expansion εv ≈ 0.97 S [2], i.e. εv ∝ C.

2 Analytical solution of the diffusion equation for constant diffusivity

First, let us consider the case of constant diffusivity, D = D0. As shown by Carslaw and Jaeger ([9], Section 2.7), the concentration profile, C(z) resulting from the boundary condition for a semi-infinite body is given by

$$\frac{C(\zeta ,\tau )}{{C_{0} }} = {\text{erfc}}\left[ {\frac{\zeta }{2}} \right] - \exp [\zeta \sqrt \tau + \tau ]{\text{erfc}}\left[ {\frac{\zeta }{2} + \sqrt \tau } \right]$$
(3)

At the surface, z = 0:

$$C(0,\,\,\tau )/C_{0} = 1 - \exp [\tau ]{\text{erfc[}}\sqrt \tau ],$$
(4)

with the normalized dimensionless time τ and normalized depth coordinate ζ, defined by

$$\tau = \frac{{h^{2} }}{{D_{0} }}t\;;\;\;\zeta = \frac{z}{{\sqrt {D_{0} \,t} }}$$
(5)

For the ratio C(ζ,τ)/C(0,τ)

$$\frac{C(\zeta ,\tau )}{C(0,\tau )} = \frac{{\text{erfc}\left[ {\frac{\zeta }{2}} \right] - \exp [\zeta \sqrt \tau + \tau ]\text{erfc}\left[ {\frac{\zeta }{2} + \sqrt \tau } \right]}}{{1 - \exp [\tau ]\text{erfc} [\sqrt \tau^{{}} ]^{{}} }}$$
(6)

two limit cases are of special interest. At very short times, τ → 0, we obtain by a series expansion with respect to τ

$$\frac{C(\zeta ,0)}{C(0,0)} = \exp \left[ { - \frac{{\zeta^{2} }}{4}} \right] - \frac{\sqrt \pi }{2}\zeta \,{\text{erfc}}\left[ {\frac{\zeta }{2}} \right]$$
(7)

At very long times τ → ∞, only the first term on the right-hand side of Eq. (3) remains finite. Consequently, we obtain the well-known solution for constant surface concentration:

$$\frac{C(\zeta ,\infty )}{C(0,\infty )} = {\text{erfc}}\left[ {\frac{\zeta }{2}} \right]$$
(8)

These limit cases are plotted in Fig. 1a. The depths at which these limit distributions decrease to C(ζ,τ)/C(0,τ) = 1/2 are

$$\begin{aligned} \zeta_{1/2} &= 0.6995_{{}} \;\quad \quad \quad \;\text{for}\quad \tau \to 0 \hfill \\ \zeta_{1/2} &= 0.9538^{{}} \;\;\quad \quad \quad \text{for}\quad \tau \to \infty \hfill \\ \end{aligned}$$
(9)
Fig. 1
figure 1

a Concentration profiles for limit cases derived from the analytical solution of diffusion, Eqs. (7, 8) with constant diffusivity, D = D0. b Water uptake according to Eq. (11). The negative sign at ΔmC stands for the decrease of the water content

Full size image

The areas under the curves define the water uptake in normalized time and depth units

$$m_{C} (\tau ) = \int\limits_{0}^{\infty } {C(\zeta ,\tau )\;d\zeta } = C(0,\tau ) \times \left\{ {\begin{array}{*{20}c} {\frac{\sqrt \pi }{{2_{{}} }}} & {{\text{for}}\quad \tau \to 0} \\ {\frac{{2^{{}} }}{\sqrt \pi }} & {{\text{for}}\quad \tau \to \infty } \\ \end{array} } \right.$$
(10)

For etching tests it is of advantage to know the amount of water ΔmC, when a layer of thickness δ has been removed from the surface. In this case it holds

$$\begin{aligned} & \Delta m_{C} (\delta ,\tau ) = - \int\limits_{0}^{\delta } {C(\zeta ,\tau )\;d\zeta } = \\ & = - m_{C} (\tau )\left\{ {\begin{array}{*{20}l} {\tfrac{1}{\sqrt \pi }\delta \exp [ - \tfrac{1}{4}_{{_{{}} }} \delta^{2} ] + \,\text{erf}[\tfrac{1}{2}\delta ]_{{}} - \tfrac{1}{2}\delta^{2} \,\text{erfc}[\tfrac{1}{2}\delta ]_{{}} } & {\text{for}\;t \to 0}\\ {(1 - \exp [ - \tfrac{1}{4}\delta^{2} ]) + \frac{\sqrt \pi }{2}^{{}} \delta \,\text{erfc}[\tfrac{1}{2}\delta ]} & {\text{for}\;t \to \infty } \\ \end{array} } \right. \\ \end{aligned}$$
(11)

The results from (11) are shown in Fig. 1b. For the thickness removal d in normal thickness unit we have to replace δ by d/√(D0t).

3 Solutions under swelling conditions

3.1 Stress enhanced diffusion

The diffusivity as a function of stress is commonly expressed by the hydrostatic stress component, σh. The diffusivity for the case of stress-enhanced diffusion is given by the following equation [10]

$$D = D_{0} \exp \left[ {\sigma_{h} \frac{{\Delta V_{w} }}{RT}} \right]$$
(12)

where D0 denotes the value of the diffusivity in the absence of a stress. T is the absolute temperature in K; ∆Vw is the activation volume for stress-enhanced diffusion and R is the universal gas constant.

The hydrostatic stress term caused by swelling stresses is

$$\sigma_{h} = - \frac{{2{\kern 1pt} \,\kappa \,E}}{9(1 - \nu )}k\,C,\quad \kappa \cong 0.97,$$
(13)

where E is Young’s modulus, ν Poisson’s ratio, and k is the equilibrium constant of the silica/water reaction given for temperatures < 500 °C by k = S/C (C = molecular water concentration, S = hydroxyl concentration).

According to Eq. (13) the swelling stress depends linearly on the water concentration, σh ∝ C. The saturation value of σh,sw for C = C0 is in the following considerations denoted as σh,0. In order to allow short expressions, the exponential term in Eq. (12) may be abbreviated by

$$\sigma_{h} \frac{{\Delta V_{w} }}{RT} = \sigma_{h,0} \frac{{\Delta V_{w} }}{RT}\frac{C}{{C_{0} }} \equiv \alpha \,C,\quad \alpha \equiv \frac{{\sigma_{h,0} }}{{C_{0} }}\frac{{\Delta V_{w} }}{RT}$$
(14)

3.2 Solution based on a perturbation set-up by Singh

By use of the Boltzmann substitution

$$\lambda = \frac{z}{2\sqrt t } \Rightarrow \lambda = \sqrt {D_{0} } \frac{\zeta }{2}$$
(15)

an ordinary differential equation results from Eq. (1)

$$2\lambda \frac{dC}{d\lambda } + \frac{d}{d\lambda }\left[ {D\frac{dC}{d\lambda }} \right] = 0 \Rightarrow \zeta \frac{dC}{d\zeta } + 2\frac{d}{d\zeta }\left[ {\frac{D}{{D_{0} }}\frac{dC}{d\zeta }} \right] = 0$$
(16)

Gardner [11] and later Singh [12] showed that this equation can be solved if the diffusion coefficient fulfills an exponential relation

$$D = \exp (\alpha C + \beta )$$
(17)

with constant coefficients α and β. This result is used in bottom mechanics [13] where the diffusivity depends on the water concentration, too.

The solution based on a perturbation ansatz reads

$$D = c_{1} \frac{{\sqrt {\pi D_{0} } }}{2}{\text{erf}}\frac{\zeta }{2} + c_{2}$$
(18)

For the swelling problem the condition (17) is fulfilled since

$$D = D_{0} {\kern 1pt} \,\exp (\Delta V_{w} \sigma_{h} /RT) = D_{0} \exp (\alpha C),$$
(19)
$$\beta = \ln D_{0} \;,\;\alpha C = \frac{{\Delta V_{w} \sigma_{h} }}{RT}$$
(20)

Combining Eqs. (18) and (19) yields

$$c_{1} \frac{{\sqrt {\pi D_{0} } }}{2}{\text{erf}}\frac{\zeta }{2} + c_{2} = D_{0} \exp (\alpha C)$$
(21)

and from this the water concentration results as a function of depth z and time t

$$C = \frac{1}{\alpha }\ln \left( {c_{1} \frac{{\sqrt {\pi /D_{0} } }}{2}{\text{erf}}\frac{\zeta }{2} + \frac{{c_{2} }}{{D_{0} }}} \right)$$
(22)

For z → ∞ it must hold C → 0, D → D0. This condition gives with erf[∞] = 1

$$c_{1} \frac{{\sqrt {\pi /D_{0} } }}{2} + \frac{{c_{2} }}{{D_{0} }} = 1\;\; \Rightarrow \quad c_{2} = D_{0} - c_{1} \frac{{\sqrt {\pi \,D_{0} } }}{2}$$
(23)

Replacing c2 in Eq. (22) results in the solution

$$C = \frac{1}{\alpha }\ln \left( {c_{1} \frac{{\sqrt {\pi /D_{0} } }}{2}{\text{erf}}\frac{\zeta }{2} + 1 - c_{1} \frac{{\sqrt {\pi /{\kern 1pt} D_{0} } }}{2}} \right)$$
(24)

or with the complementary error function erfc(x) = 1 − erf(x):

$$C = \frac{1}{\alpha }\ln \left( {1 - c_{1} \frac{{\sqrt {\pi /D_{0} } }}{2}{\text{erfc}}\frac{\zeta }{2}} \right)$$
(25)

3.2.1 Increasing surface concentration (mass-transfer condition)

A solution for the surface conditions by Eq. (2) cannot result from Singh’s procedure. This can even be seen from the application of the Boltzmann substitution. In terms of the normalized time and depth coordinate by (5), the substitution λ in Eq. (15) is only dependent on the depth coordinate ζ and not the time τ. Consequently the applicability to time-dependent diffusion effects given by Eq. (2) is strongly restricted. Nevertheless, this solution is appropriate for treating the limit case τ → ∞, i.e. for the condition of fixed surface concentrations.

3.2.2 Constant surface concentration

For very long diffusion times, the surface water concentration, C(0), tends asymptotically to the saturation value C0. In order to compute the limit case for t → ∞. Specimen soaked in water vapour for very long times are assumed to show constant surface water concentration C(z = 0) = C0. In this case we obtain from (25)

$$C_{0} = \frac{1}{\alpha }\ln \left( {1 - c_{1} \frac{{\sqrt {\pi /D_{0} } }}{2}} \right)$$
(26)

the constant c1 as

$$c_{1} = 2\sqrt {D_{0} /\pi } (1 - \exp [\alpha C_{0} ])$$
(27)

The result for the concentration is then

$$\frac{C}{{C_{0} }} = \frac{1}{{\alpha C_{0} }}\ln \left( {1 - (1 - \exp [\alpha C_{0} ])\,\text{erfc}\left[ {\frac{\zeta }{\text{2}}} \right]} \right)$$
(28)

Water profiles computed via Eq. (28) are shown in Fig. 2a for different parameters αC0. Figure 2b shows a comparison of the analytical solution Eq. (28) as the black curve and the numerical results according to [2] as the red curve, both for αC0 = − 3. The small differences may be the consequence of the finite depth interval of 5 × δ1/2 (i.e. the concentration boundary condition at this location prescribed by C(ζ = 5 × δ1/2) = 0) that had to be used in the numerical program NDSolve by Mathematica [14].

Fig. 2
figure 2

a Effect of swelling on the water profiles, b comparison of Eq. (28) with numerical solution from [2], given by the black and red curve, respectively

Full size image

The depth ζ1/2 at which the distributions of Fig. 2a decreases to C(ζ)/C0 = 1/2 is

$$\zeta_{1/2} = 2\,\text{erf}^{ - 1} \,\left[ {\infty , - \frac{1}{{1 + \exp [\tfrac{1}{2}\alpha {\kern 1pt} C_{0} ]}}} \right]$$
(29)

Numerical values are compiled in the second column of Table 1.

Table 1 Data of water profiles and water uptake obtained for saturation conditions, τ → ∞
Full size table

Finally, we determined water uptake by integrating the swelling profiles of Fig. 2a numerically with the result

$$m_{C} (\tau ) = \int\limits_{0}^{\infty } {C(\zeta ,\tau )\;d\zeta }$$
(30)

with results given in the third column of Table 1 for τ → ∞. The value for αC0 = 0 is 2/√π ≅ 1.128 as given by Eq. (10).

The decrease of the water by surface removal d is shown in Fig. 3 in terms of the normalized quantity δ

$$\delta = \frac{d}{{\sqrt {D_{0} \,t} }}$$
(31)
Fig. 3
figure 3

Change of water uptake with normalized surface removal depth δ as a function of swelling parameter αC0 normalized on the uptake for αC0 = 0

Full size image

The depths at wich half of the water content is removed, δ1/2, is given in the fourth column of Table 1 for a few values of αC0.

4 Discussion

Results of computations are shown in Fig. 4a, b. The bold curve in Fig. 4a a is the limit case for saturation conditions for τ → ∞, given by Eq. (25) and the thin curves represent the numerical solution of Eqs. (1) and (2). The widths of the profiles allow an effective diffusivity to be determined according to Davis and Tomozawa [15]

$$C = C_{0} \,{\text{erfc}}\left( {\frac{z}{{2\sqrt {D_{eff} \,t} }}} \right)$$
(32)
Fig. 4
figure 4

a Diffusion profiles as a function of time without and with consideration of swelling stress effect on diffusivity (thin curves computed numerically, bold curves: analytical solutions), b effective diffusivity Deff vs. normalized time τ, c effective diffusivities measured by Oehler and Tomozawa [16] at 250 °C and 39 bar water vapour pressure

Full size image

In [15] the effective diffusivity, Deff, is determined from the width at C/C0 = 1/2 interpreting it as the diffusion lengths \(\;\sqrt {D_{eff} \,t}\). This effective diffusivity Deff is usually determined in investigations on water diffusion in silica [15,16,17,18].

Measurements of effective diffusivities by Oehler and Tomozawa [16] at 250 °C and 39 bar water vapour pressure are introduced by circles, see Fig. 4c. The curve is introduced as a guideline for the eyes. A clear decrease with time is evident. This effect was discussed in [2] in terms of swelling stresses.

5 Summary

Water present at silica surfaces reacts with the silica network under generation of hydroxyl that causes a volume expansion of the thin water-affected zone. Since free expansion is restricted by the bulk material, compressive swelling stresses result which yield in a reduced diffusivity. Stress affected diffusion was considered and resulting saturation profiles were derived by the analytical procedure of Gardner [11] and Singh [12]. The calculations show that the compressive swelling stresses result in steeper profiles caused by a decreased diffusivity compared to stress-free state. In order to show the occurrence of reduced diffusivity, we compared the predicted water profiles with measurements by Öhler and Tomozawa [16]. The predictions made by computations show the same tendency as experimental results from literature which couldn’t be interpreted sufficiently so far.