*2.2. Theoretical Analysis*

This paper analyzes the hydrogen venting in the Sao Francisco basin analytically. This reveals the fundamental controls on the pulsing venting, but numerical simulations will ultimately be needed to gain a full understanding of the system. Our analysis is thus preliminary and preparatory. All symbols are defined in Table 1. Parameter values for the problem at hand are also given in the table, along with the method of their calculation or references for the values selected. Relationships given in the table indicate important physical dependencies. We consider both temperature and pressure wave transmission into the subsurface but concentrate on pressure transmission because we find that daily temperature changes penetrate only a few 10 s of centimeters and therefore cannot modulate hydrogen concentrations to a depth of one meter. Our method is to calculate the change in gas volume from either thermal expansion or compression and integrate this volume change to determine vertical (1D) gas velocity.

### Thermal and Pressure Wave Propagation into the Subsurface

Surface temperature and pressure changes diffuse into the subsurface at rates and depths controlled by the subsurface thermal and hydraulic diffusivity, respectively. The diffusion equations both combine the heat and pressure flux equations with conservation of heat and mass and are identical in form and most clearly described in parallel. These equations describe how cyclic variations of temperature and pressure propagate from the surface at z = 0 into the subsurface (z negative):

The flux of heat described by Fourier's law and the flux of gas by Darcy's law are:

$$j = -K \frac{\partial T}{\partial z} 2^{\frac{1}{2}} \tag{1}$$

$$V = -\frac{k'}{\mu} \frac{\partial p'}{\partial z} \mathbf{\hat{z}} \tag{2}$$

Conservation of mass requires:

$$\begin{array}{l}\rho\_T c\_T \frac{\partial T}{\partial t} = -\nabla \bullet j = K\_T \frac{\partial^2 T}{\partial z^2} \\ \frac{\partial T}{\partial t} = \kappa\_T \frac{\partial^2 T}{\partial z^2} \\ \kappa\_T = \frac{K\_T}{\rho\_T c\_T} \end{array} \tag{3}$$

$$\begin{cases} \beta\_m + \phi \beta\_f \Big| \frac{\partial p'}{\partial t} = -\nabla \bullet V = \frac{k'}{\mu\_{air}} \frac{\partial^2 p'}{\partial z^2} \\ \frac{\partial p'}{\partial t} = \kappa\_p \frac{\partial^2 p'}{\partial z^2} \\ \kappa\_p = \frac{k'}{\mu\_{air} (\beta\_m + \phi \beta\_{air})} \cong \frac{k'}{\mu\_{air} \phi \beta\_{air}} \end{cases} \tag{4}$$

In (4) we have used the fact that the compressibility of air is much greater than the compressibility of the soil matrix.

*Geosciences* **2020**, *10*, 149

The final temperature (3) and pressure (4) diffusion equations differ only in their diffusivity parameters κ. We seek the subsurface solution for *T*(*<sup>z</sup>*,*<sup>t</sup>*) and *p*'(*<sup>z</sup>*,*<sup>t</sup>*) for *z* = 0 to −∞. Let α represent either *T* or *p'* and αo the amplitude variation imposed at the surface. Defining *t* = *t*/*P*, where P is the period of the harmonic variation imposed at the surface, the solution to Equations (3) and (4), (Carslaw and Jager, 1959. p.65 [24]) is:

$$\alpha(z,\tilde{t}) = a\_o e^{z/\delta} \cos(2\pi \tilde{t} + z/\delta) \tag{5}$$

where

$$
\delta\_{\alpha} = \sqrt{\frac{\kappa\_{\alpha} P}{\pi}} \tag{6}
$$

Because we want the rate of venting, we take the time derivative of (5):

$$\frac{\partial \alpha}{\partial t} = \frac{1}{P} \frac{\partial \alpha}{\partial \tilde{t}} = -\frac{2\pi}{P} \alpha\_o e^{z/\delta\_a} \sin \left(2\pi \tilde{t} + z/\delta\_a\right) \tag{7}$$

Defining *z* = *<sup>z</sup>*/δα, we can integrate (5) from some depth *zb* to the surface (or sensor depth) to obtain the total rate of volume change from grea<sup>t</sup> depth to the subsurface (or up to sensor depth). By conservation of mass this must equal the gas flux, *Vz*=<sup>0</sup>(*t*), out the top surface (or past the sensors). The sign depends on whether the forcing is by pressure or temperature change. The gas fluxes are:

$$V\_{\Xi=0}^{p'}(t) = -\frac{2\pi}{P} \alpha\_{p'o} \phi \beta\_{p'} \delta\_{p'} \int\_{\Xi'=\Xi\_b}^{0} \varepsilon^{\overline{z}'} \sin(2\pi \overline{t} + \overline{z}') d\overline{z}' \tag{8}$$

$$V\_{\Xi=0}^{T}(t) = \frac{2\pi}{P} \alpha\_{To} \phi \beta\_T \delta\_T \int\_{\substack{\Xi'=\Xi\_b}}^{0} e^{\overline{\Xi'}} \sin(2\pi \overline{t} + \overline{z}') d\overline{z}' \tag{9}$$

Here α*o* has been replaced by <sup>α</sup>*p'o* or α*To*, and the sign is such that there is volume expansion if the change in pressure is negative or the change in temperature is positive. Note the prime on *z* simply indicates that it is an integration variable.

If the depth of integration, *zb*, is sufficiently large that the underlying pressure or temperature variations are negligible, for example *zb* = −10, the normalized integral, 0 *<sup>z</sup>*<sup>=</sup>−10 *ez* sin(<sup>2</sup>π*<sup>t</sup>* + *z*)*dz* varies from −0.709 to +0.709 as *t* varies over its sinusoidal cycle. The phase shift between the temperature or

pressure forcing applied at the surface (*z* = 0) and the venting flux *V* are both maximum for this deep integration. If *zb* is shallow and the full subsurface thermal or pressure wave is not captured in the integration, the phase shift between the pressure or temperature and the venting flux *V* at the *z* = 0 is less and the rate of venting is less.

The maximum and minimum values of the harmonic gas flux at the surface are:

$$V\_a^{\text{max}/\text{min}}(z=0) = \pm \frac{(0.709)2\pi}{P} a\_o \phi \beta\_a \delta\_a \tag{10}$$

Since the gas flux is linear, we can add the deep *<sup>H</sup>*2-bearing gas flux, *Vdeep*, to (6) to obtain the total gas flux, *Vtotal <sup>Z</sup>*=<sup>0</sup>(*t*):

$$V\_{z=0}^{\text{total}}(t) = V\_{\text{deep}} - \frac{2\pi}{P} \alpha\_0 \phi \beta\_0 \delta \int\_{\substack{\tilde{z}' = \tilde{z} \\ \tilde{z}' = -10}}^{\tilde{z}} \tilde{e}' \sin(2\pi \tilde{t} + \tilde{z}') d\tilde{z}' \tag{11}$$

The H2 emissions could be pulsating if:

$$V\_{dexp} < \frac{2\pi}{P} \alpha\_o \phi \beta\_a \delta\_a(0.709) \tag{12}$$


**Table 1.** Parameter values in equations derived above.


**Table 1.** *Cont.*

If condition (12) is satisfied, the T-driven or p'-driven flow into the ground will exceed the base leakage rate of deep *<sup>H</sup>*2-bearing gas (*Vdeep*), *H*2-free air could periodically surround the sensors, and the sensors could then detect low *H*2 concentrations. Since *Vdeep* = 0.05 to 1.2 m<sup>3</sup> m<sup>−</sup><sup>2</sup> d−1, the question is whether the right hand side of (12) can exceed this value.
