What Pulsating H₂ Emissions Suggest about the H₂ Resource in the Sao Francisco Basin of Brazil

Abstract

Proterozoic sedimentary basins very often emit natural hydrogen gas that may be a valuable part of a non-carbon energy infrastructure. Vents in the Sao Francisco Basin in Brazil release hydrogen to the atmosphere mainly during the daylight half of the day. Daily temperature and the regular daily tidal atmospheric pressure variations have been suggested as possible causes of the pulsing of H₂ venting. Here, we analyze a ~550 m-diameter depression that is barren of vegetation and venting hydrogen mainly at its periphery. We show that daily temperature changes propagated only ~1/2 m into the subsurface and are thus too shallow to explain the H₂ variations measured at 1-m depth. Pressure changes could propagate deeply enough, and at the depth at which the cyclic variations are measured hydrogen concentration will have the observed phase relationship to atmospheric pressure changes provided: (1) the pressure wave is terminated by geologic barriers at about 25% of its full potential penetration distance, and (2) the volume of gas in the vents is very small compared to the volume of gas tapped by the venting. These constraints suggest that there is a shallow gas reservoir above the water table under the ~550 m-diameter barren-of-vegetation depression. The 1D-analytical and finite-element calculations presented in this paper help define the hydrogen system and suggest the further steps needed to characterize its volume, hydrogen flux and resource potential.

1. Introduction

Hydrogen gas may become an important part of the zero-carbon energy economy because its combustion produces only water vapor and no CO₂. Hydrogen is generally considered a vector of energy because it is an agent of energy transfer that needs to be manufactured from some other material, such as hydrocarbons or water, and thus carries energy from other sources. Several studies have shown, however, that there are natural sources of hydrogen that could be an important resource (Ward, 1933 [1]; Goebel et al., 1984 [2]; Newell et al., 2007 [3]; Sherwood-Lollar et al., 2014 [4]; Prinzhofer and Deville, 2015 [5]; Guélard et al., 2017 [6]; Prinzhofer et al., 2018 [7]). The Earth is naturally expelling native hydrogen through still-poorly determined physico-chemical processes.

As at the start of the petroleum era, we currently know only that hydrogen is seeping out of rocks and soils in many locations. Studies so far have focused on the chemistry of the seeps that have been discovered worldwide (Larin et al., 2014 [8]; Zgonnik et al., 2015 [9]; Deville and Prinzhofer, 2016 [10]; Prinzhofer et al., 2019 [11]) and the origin of the hydrogen (Larin, 1993 [12]; Gilevska, T., 2007 [13]; Milesi et al., 2015 [14]; Vacquand et al., 2018 [15]; Truche et al., 2018 [16]). The many seeps that have been found indicate the operation of an active hydrogen system, and transient accumulations of hydrogen at relatively shallow depth are known (Mali, Prinzhofer et al., 2018 [7], USA, Goebel et al., 1984 [2]; Newell et al., 2007 [3]; Guélard et al., 2017 [5]). But we know little about the scale and economic potential of the hydrogen systems that are involved.

The most important parameter for gauging the magnitude of the hydrogen resource is the magnitude of the hydrogen (H₂) emission because it immediately suggests the probability of hydrogen accumulation under subsurface structures, and also indicates the sustainability of this source of energy. To better constrain the magnitude of the hydrogen emission, a permanent monitoring station has been installed at 16 °S in the dry, hot Sao Francisco Basin of Brazil. Venting is being monitored there with sensors that measure H₂ concentrations at ~1 m depth around the perimeter of a 550 m diameter vent (Prinzhofer et al., 2019) [11]. The data indicate daily pulses of hydrogen emission. The maximum concentration of H₂ in the venting gases usually occurs in the middle of the day, although some sensors have hydrogen concentration peaks at the beginning of the day, or even during the night. It appears that these variations are linked to complex external parameters such as biologic activity in the soils, temperature, atmospheric pressure, Earth tides, etc.

Importantly, the correlation between H₂ concentrations and atmospheric pressure is not linear, but presents a kind of hysteresis in which there is a delay between the pressure variation in free air, the hydrogen concentration measured by the sensors at 1 m depth, and the gas movements induced in the porous soil. The soil generally has very low water saturation above the water table several meters below the surface. This paper investigates whether the variation of hydrogen concentrations in soils could be caused by diurnal/nocturnal variations in temperature or atmospheric pressure.

H₂ gas is venting in the Sao Francisco Basin from slight topographic depressions that are circular, barren of vegetation, and aligned along a major fault. The water table is at 3 to 5 meters depth. The hydrogen venting from the 600 × 500 m study area is not continuous, but pulses on a daily cycle wherein H₂ is detected by shallow (1 m depth) sensors for only about half the day, usually centered on high noon. Prinzhofer et al. (2019) [11] have estimated the emission rate is between 7000 and 178,000 m³ H₂ per day, based on the venting at Sensor 9 (Figure 1) where the volume of H₂ venting is greatest and the concentration of H₂ in the venting gas is sometimes above 1000 ppm. Hydrogen is thought to be generated in basins from the reduction of water attending ferrous to ferric iron oxidation (Welhan and Craig, 1979 [17]; Neal and Stranger, 1983 [18]; Abrajano et al., 1990 [19]; Sano et al., 1993 [20]; Charlou et al., 2002 [21]) and/or NH₄⁺ decomposition at depths where the temperature exceeds 200 to 250 °C (Guélard et al., 2018) [22]. The H₂ concentration was found to be 50% to 80% (with the rest N₂ and several percent CH₄) in a deep drill hole in the Sao Francisco basin, and measurements at other sites in the basin found a source H₂ concentration ranging from 50% to ~100%. This suggests that the undiluted hydrogen source at the vent we analyze has a hydrogen concentration of at least 50%.

We investigate whether, and under what circumstances, daily changes in either surface temperature or atmospheric pressure can cause gas flow to reverse (e.g., air to flow into the ground) for half the day and thus potentially explain why H₂ is detected for only half a day. The initial method of analysis is to compute, for both temperature and pressure changes, the change in gas density as a function of depth, the variations in the rate of gas flow produced by these density changes, and the conditions under which the density-change-driven-flow could cause the flow of H₂-free air into the vent subsurface for about half the day. It is found that variations in atmospheric pressure could account for the diurnal variations in H₂ detection, but temperature variations cannot. The pressure-driven changes in hydrogen concentration at 1 m depth are then investigated quantitatively by comparing the timing and magnitude (phase relations between) of the calculated and observed modulation in hydrogen concentration. It is found that pressure variations must be blocked by a barrier of some kind when they have propagated about 25% of the distance they could potentially propagate into a gas-filled reservoir, and that the gas volume of the reservoir tapped needs to be at least 1000 times greater than the gas volume in the vent. These insights, obtained by simple 1D analytical and finite element analyses, suggest how the Sao Francisco hydrogen system is operating, how the system could be tested by future drilling, and the kind of future 3D modeling that will be required to define the operation of this hydrogen system.

2. Analysis

2.1. Site Characterisics

2.1.1. The Study Site

As illustrated in Figure 1, H₂ is venting from a 600 × 500 m depression that is barren of vegetation. One possibility is that the reduced conditions caused by hydrogen venting have killed the vegetation; another is that periodic flooding of the small depressions is he cause. To be neutral, we label the area A_barren. The study area is the large barren area shown in Figure 1 with A_barren = 300,000 m². It is venting hydrogen at an estimated rate of between 7000 and 178,000 m³ H₂ per day (Prinzhofer et al., 2019) [11]. Assuming the concentration of the source is 50% H₂, the venting rate at depth, Q_H2, is between 14,000 and 356,000 m³ per day. If the flow is uniform across the vent area, the gas flux is V_deep = Q_H2/A_barren = 0.05 to 1.2 m d⁻¹. The average concentration of H₂ in the discharge measured at ~1 m depth is ~100 ppm (Prinzhofer et al., 2019 their Figure 2) [11].

2.1.2. Daily Variations in T

The temperature ranges from 15 to 35 °C in the afternoons. Because the bare surface will likely be heated above the atmospheric temperature, the range of ground surface temperature will likely be greater than the average atmospheric temperature range.

2.1.3. Daily Variations in Atmospheric Pressure

The daily variations in atmospheric pressure are quite regular at the site and are available from a local meteorological station. Pressure peaks at about 10:00 local solar time, and the steepest part of the subsequent pressure decrease occurs at about 13:00 local solar time. The pressure change over the day is between 4 and 8 mbar. Figure 2A compares the atmospheric pressure variations on 14 August 2018 to the [H₂] concentrations measured at 1 m depth by Sensor 6 on the same day. Peak [H₂]_sensor-6 occurs at ~13:00, coinciding with the maximum in the rate of atmospheric pressure decline.

Figure 2B separates out the diurnal and semidiurnal sinusoidal components of the atmospheric pressure change amplitudes of 1.5 and 1.2 mb and phase lags of 16 and 14 h, respectively. When summed these components replicate the observed pressure variations very closely (compare yellow and blue curves in Figure 2A).

Figure 2C shows the phase relationships between atmospheric pressure and the hydrogen concentration measured at Sensor 6. As in Figure 2A, hydrogen is detected when the atmospheric pressure is decreasing, but the figure shows a clear hysteresis between pressure change and the hydrogen measurements that is diagnostic of pressure wave transmission into the subsurface. By hysteresis we mean that hydrogen concentrations measured in Sensor 6 when the atmospheric pressure is increasing are different from those measured when the atmospheric pressure is decreasing. In this case the contrast is dramatic since the hydrogen concentration is zero when pressure is increasing and positive when it is decreasing. The red circle in Figure 2C emphasizes that the maximum [H₂] is measured at 13:00 when the surface pressure is decreasing at its maximum rate (Figure 2A).

The regular variations in atmospheric pressure are caused by atmospheric tides. At low latitudes, away from the much larger pressure variations associated with shifts of the polar front jets, the atmospheric tidal variations in pressure are typically mainly semi-diurnal with amplitudes of ~1.6 mb (Le Blancq, 2011 [23]). At our site the tidal atmospheric changes (2 mb amplitude) are similar in magnitude but have a strong diurnal component. We do not know why the daily atmospheric pressure variation at our site are more strongly diurnal than the norm.

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 in equations derived above.

Param.	Units	Definition	Value		Method/Comment
αoT	K	Amplitude harmonic Tz = 0(t)	30		Ground changes > atmospheric.
αop’	Pa	Diurnal amplitude p’z = 0(t) (lag 16 h) Semidiurnal p’z = 0(t) (14 h phase lag)	150 120		Atmospheric tide has two harmonic components (Figure 2).
βT	-	Thermal expansion of air	3.3 × 10−3		=1/T[K] = 1/298 for 25 °C.
βp’	-	Volumetric compression of ideal gas	10−5		=1/p = 1/105 at atm. Pressure.
b	m	Depth of system Equation (17)	~20 m		[H2] changes only near the surface.
bbox	m	Depth extent box model Equation (13)	>500 m		Large bbox needed for sufficient Vair
cs	J kg−1 K−1	Heat capacity of solid matrix	837
cair	J kg−1 K−1	Heat capacity of water-saturated air	4380
δT	m	Thermal skin depth	ϕ	δT	$=\sqrt{\frac{{\kappa }_{T}P}{\pi }},P=3600\times 24s$
			0.1	0.107
			0.2	0.102
			0.4	0.098
${\delta }_{p’}^{1d}$	m	Pressure skin depth for 1 day period	ϕ	δp’	$=\sqrt{\frac{{\kappa }_{p’}P}{\pi }},P=3600\times 24s$. Reduced by $\sqrt{2}$ for semidiurnal.
			0.1	56.1
			0.2	39.7
			0.4	28.0
D	m2 s−1	Diffusion constant of H2 in large excess of air	0.756 × 10−4 6.45 m2 d−1		By comparison DCO2= 0.208 × 10−4. Values from web.
DE	m2 s−1	Effective diffusion constant in vent	ϕ	DE	$=\frac{D\varphi }{\tau }$ values assume τ = 2. Table units are m2 per day.
			0.1	0.323
			0.2	0.645
			0.4	1.29
ϕ	-	Porosity of vent subsurface	~0.4
[H2]	-	Volume fraction of H2 in air
[H2]1m	-	Volume fraction of H2 in air at 1 m depth		<102–103 ppm
[H2]deep	-	Fraction H2 deep in the vent	0.5 to 1		Deep H2 concentration ~1 from diverse observations (see text).
$\left[{\overline{H}}_2\right]$	-	H2 as fraction of [H2]deep			=[H2]/[H2]deep.
j	J m−2 s−1	Heat flux	calculated
k’	m2	Permeability (10−12 m2 = 1 Darcy)	10−12		Near surface sandy soil is assumed to be very permeable.
KT	W m−1 K−1	Thermal conductivity of media	ϕ	KT	Calculated at 25 °C with fabric theory mix of Ks and Kair with fabric mixing parameter = 0.3, see Luo et al. (1994) [25].
			0.1	0.843
			0.2	0.685
			0.4	0.498
Ks	W m−1 K−1	Matrix thermal conductivity	4		Sandy soil with some clay at 25 °C.
Kair	W m−1 K−1	Air thermal conductivity	0.0261		Air at 25 °C. T Dependence weak.
κT	m2 s−1	Thermal dispersivity of media	ϕ	κT	$=\frac{K_T}{{\rho }_T{c}_T}$.
			0.1	4.2 × 10−7
			0.2	3.8 × 10−7
			0.4	3.6 × 10−7
κp’	m2 s−1	Hydrodynamic dispersivity of H2-filled media for 1 Darcy vent permeability	ϕ	κp’	$\cong \frac{k^{\prime }}{{\mu }_{H_2}\varphi {\beta }_{p’}}$,κp’ for air is half H2.
			0.1	0.114
			0.2	0.057
			0.4	0.029
Npe	-	Peclet nbr (Equation (17))	>7.72		$=\frac{Vb}{D_E}=\frac{(0.05)(100)}{0.648}=7.72$
μair, μH2	kg m−1 s−1	dynamic viscosity of air	μair = 1.78 × 10−5 μH2 = 8.74 × 10−6		μair = νair ρair.
νair	m2 s−1	kinematic viscosity of air	1.48 × 10−5		Regression: νair = −1.1555(10−14)T3 + 9.5728(10−11)T2 + 3.7604(10−8)T − 3.4484(10−6), from WWW.
P	s	Period of P harmonic	24 × 3600		T surface change, tidal p’ change.
$\overline{P}$	-	Period normalized to 1 day	1		P/(1 day).
p	Pa	Atmospheric pressure	105
P’	kg m s−2	Perturbation from atmos. pressure
p’o	kg m s−2	Amplitude atm pressure change	200 to 400		Changes due to atmospheric tide.
ρT	kg m−3	Density of media	ϕ	ρT	=ρs(1 − ϕ) + ϕρair.
			0.1	2430
			0.2	2160
			0.4	1620
ρs	kg m−3	Density of solid matrix	2700
ρair	kg m−3	Density of air	1.1		Density of air at 25 °C.
ρTcT	J kg−1 K−1	Volumetric heat capacity air-saturated soil	ϕ	cT	=ρscs(1 − ϕ) + ϕρaircair.
			0.1	2 × 106
			0.2	1.8 × 106
			0.4	1.4 × 106
T	K	Temperature in degrees K	~298
To	K or °C	Amplitude of imposed surface T(t)	>20		Can be greater if surface T > T air.
τ	-	Tortuosity of sediment pores	2
τD	d	Time normalization in Equation (16)			$=\frac{b^2\phi }{D_E}$.
Vair	m3 m−2 d−1	Flux of air containing H2			Vair = Vcyclic + Vdeep.
Vcyclic	m3 m−2 d−1	Flux of air driven only by changes p
Vdeep	m3 m−2 d−1	Deep (source) H2 gas flux assuming 100% H2 concentration	0.05 to 1.2		Estimated by Prinzhofer (2019) [11], assuming uniform venting in barren zone.
VTmax	m d−1	Maximum flux into subsurface from T(t)	ϕ	VTmax	$=\frac{2\pi }{(P=1d)}{\alpha }_{oT}\varphi {\beta }_T{\delta }_T(0.709)$
			0.1	0.005
			0.2	0.009
			0.4	0.018
Vp’max diurnal	m d−1	Maximum flux into subsurface from diurnal p’(t), αop’ = 150 Pa. From Equation (6)a. A 1 darcy subsurface permeability is assumed.	ϕ	Vp’max	$=\frac{2\pi }{(P=1d)}{\alpha }_{op’}\varphi {\beta }_{p’}{\delta }_{p’}^{1d}\left(0.709\right)$
			0.1	0.037
			0.2	0.053
			0.4	0.075
Vp’max Semi-diurnal	m d−1	Maximum flux into subsurface from semidiurnal p’(t), αop’ = 120 Pa. From Equation (8). A 1 darcy vent permeability is assumed.	ϕ	Vp’max	$=\frac{2\pi }{(P=0.5d)}{\alpha }_{op’}\phi {\beta }_{p’}\frac{{\delta }_{p’}^{1d}}{\sqrt{2}}(0.709)$
			0.1	0.042
			0.2	0.060
			0.4	0.084

. 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}\widehat{z}$$

(1)

$$V=-\frac{k^{\prime }}{\mu }\frac{\partial p^{\prime }}{\partial z}\widehat{z}$$

(2)

Conservation of mass requires:

$$\begin{array}{l}{\rho }_T{c}_T\frac{\partial T}{\partial t}=-\nabla •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}$$

(3)

$$\begin{array}{l}\left({\beta }_m+\varphi {\beta }_f\right)\frac{\partial p^{\prime }}{\partial t}=-\nabla •V=\frac{k^{\prime }}{{\mu }_{air}}\frac{{\partial }^2{p}^{\prime }}{\partial z^2}\\ \frac{\partial p^{\prime }}{\partial t}={\kappa }_p\frac{{\partial }^2{p}^{\prime }}{\partial z^2}\\ {\kappa }_p=\frac{k^{\prime }}{{\mu }_{air}\left({\beta }_m+\varphi {\beta }_{air}\right)}\cong \frac{k^{\prime }}{{\mu }_{air}\varphi {\beta }_{air}}\end{array}$$

(4)

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

The final temperature (3) and pressure (4) diffusion equations differ only in their diffusivity parameters κ. We seek the subsurface solution for T(z,t) and p’(z,t) for z = 0 to −∞. Let α represent either T or p’ and α_o the amplitude variation imposed at the surface. Defining $\overline{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. p65 [24]) is:

$$\alpha (z,\overline{t})={\alpha }_o{e}^{z/\delta }\cos \left(2\pi \overline{t}+z/\delta \right)$$

(5)

where

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

(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 \overline{t}}=-\frac{2\pi }{P}{\alpha }_o{e}^{z/{\delta }_{\alpha }}\sin \left(2\pi \overline{t}+z/{\delta }_{\alpha }\right)$$

(7)

Defining $\overline{z}=z/{\delta }_{\alpha }$, we can integrate (5) from some depth ${\overline{z}}_b$ to the surface (or sensor depth) to obtain the total rate of volume change from great depth to the subsurface (or up to sensor depth). By conservation of mass this must equal the gas flux, $V_{z=0}(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_{\overline{z}=0}^{p’}(t)=-\frac{2\pi }{P}{\alpha }_{p’o}\varphi {\beta }_{p’}{\delta }_{p’}{\displaystyle \underset{\overline{z}’={\overline{z}}_b}{\overset{0}{\int }}{e}^{\overline{z}’}\sin \left(2\pi \overline{t}+\overline{z}’\right)}d\overline{z}’$$

(8)

$$V_{\overline{z}=0}^T(t)=\frac{2\pi }{P}{\alpha }_{To}\varphi {\beta }_T{\delta }_T{\displaystyle \underset{\overline{z}’={\overline{z}}_b}{\overset{0}{\int }}{e}^{\overline{z}’}\sin \left(2\pi \overline{t}+\overline{z}’\right)}d\overline{z}’$$

(9)

Here α_o has been replaced by α_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 $\overline{z}$ simply indicates that it is an integration variable.

If the depth of integration, ${\overline{z}}_b$, is sufficiently large that the underlying pressure or temperature variations are negligible, for example ${\overline{z}}_b$ = −10, the normalized integral, ${\displaystyle \underset{\overline{z}=-10}{\overset{0}{\int }}{e}^{\overline{z}}}\sin (2\pi \overline{t}+\overline{z}’)d\overline{z}’$ varies from −0.709 to +0.709 as $\overline{t}$ varies over its sinusoidal cycle. The phase shift between the temperature or pressure forcing applied at the surface ($\overline{z}$ = 0) and the venting flux V are both maximum for this deep integration. If ${\overline{z}}_b$ 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 $\overline{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_{\alpha }^{\max /\min }\left(z=0\right)=\pm \frac{(0.709)2\pi }{P}{\alpha }_o\varphi {\beta }_{\alpha }{\delta }_{\alpha }$$

(10)

Since the gas flux is linear, we can add the deep H₂-bearing gas flux, V_deep, to (6) to obtain the total gas flux, $V_{Z=0}^{total}(t)$:

$$V_{z=0}^{total}(t)=V_{deep}-\frac{2\pi }{P}{\alpha }_o\varphi {\beta }_{\alpha }\delta {\displaystyle \underset{\overline{z}’=-10}{\overset{\overline{z}}{\int }}{e}^{\overline{z}’}\sin \left(2\pi \overline{t}+\overline{z}’\right)}d\overline{z}’$$

(11)

The H₂ emissions could be pulsating if:

$$V_{deep}<\frac{2\pi }{P}{\alpha }_o\varphi {\beta }_{\alpha }{\delta }_{\alpha }(0.709)$$

(12)

If condition (12) is satisfied, the T-driven or p’-driven flow into the ground will exceed the base leakage rate of deep H₂-bearing gas (V_deep), H₂-free air could periodically surround the sensors, and the sensors could then detect low H₂ concentrations. Since V_deep = 0.05 to 1.2 m³ m⁻² d⁻¹, the question is whether the right hand side of (12) can exceed this value.

• The box model

It is instructive to consider a simple (but unrealistic) “box” model where, as atmospheric pressure changes, gas pressure changes instantly and uniformly within the entire box. If the subsurface where hydrogen is venting is permeable enough to some depth b_box, the atmospheric pressure change can be transmitted immediately, and the venting flux will be coincident with the rate of atmospheric pressure change:

$$V_{box}(\overline{z}=0)=-\varphi b_{box}\frac{\partial p}{\partial t}{\beta }_{p’}$$

(13)

For φ = 0.2, b_box = 1000 m, ${\beta }_{p’}$ = 10⁻⁵ Pa⁻¹, and −$\frac{\partial p}{\partial t}$ = 25 mb d⁻¹ = 2500 Pa d⁻¹, the gas flux from the box V_box = 5 m d⁻¹. For later discussion it is important to emphasize that the assumption of very rapid pressure transmission means that there is no phase shift between surface pressure change and surface gas flux.

• Advection-diffusion equation for [H2]

It will take some time for gas moving into the subsurface to reach the sensors at 1 m depth, and on the way, there will be mixing with the gas in the subsurface. The advection–diffusion equation describes these phenomena. The 1D flux, j_H2, of H₂ (expressed as the volume fraction of H₂ in air) by diffusion and advection is:

$$j_{H2}=-D_E\frac{\partial [H2]}{\partial z}+V_{air}[H_2]$$

(14)

where [H₂] is the volume fraction of H₂ in air. By mass conservation:

$$\frac{\partial \phi \left[H_2\right]}{\partial t}=-\nabla •j_{H2}=\frac{\partial }{\partial z}{D}_E\frac{\partial \left[H_2\right]}{\partial z}-V_{air}\frac{\partial \left[H_2\right]}{\partial z}$$

(15)

where D_E is the effective diffusion constant of H₂ in the vent, and V_air is the vertical flux of air (m³ of air passing through a plan area or 1 m² per second) in the portion of the vent considered (in this case the uppermost part). The effective diffusion constant $D_E=\frac{D\varphi }{\tau }$, where D is the diffusion constant of H₂ in air, φ is the porosity of the vent sediments, and τ is the tortuosity of the pores in the sediments, which we take to equal 2. If we make t, z and [H₂] non-dimensional by defining:

$$\begin{array}{l}t=\overline{t}\frac{b^2\phi }{D_E}=\overline{t}{\tau }_D\\ z=\overline{z}b\\ [H_2]=[{\overline{H}}_2]{[H_2]}_{deep}\\ N_{pe}=\frac{Vb}{D_E}\end{array}$$

(16)

Equation (15) becomes:

$$\frac{\partial \left[{\overline{H}}_2\right]}{\partial \overline{t}}=\frac{{\partial }^2\left[{\overline{H}}_2\right]}{\partial {\overline{z}}^2}-N_{pe}\frac{\partial \left[{\overline{H}}_2\right]}{\partial \overline{z}}$$

(17)

where N_pe is the Peclet number (ratio of advection to diffusion).

2.3. Site Analysis

2.3.1. The Minimum Criterion for [H₂]_Station-6 Modulation

Whether the pulsing H₂ leakage criterion (9) is met depends on the parameter values. The last 3 lines in the Table 1 show the calculated cyclic gas inflow rates for temperature and pressure variations at the Sao Fancisco site. The maximum temperature-driven gas influx is an order of magnitude less than the pressure-driven influx, and too small to reverse even Prinzhofer’s lowest venting estimate. The pressure-driven influx, on the other hand, could reverse the [H₂] venting if the venting were at the low end of Prinzhofer’s estimated range. For convenience we show the comparison in

Table 2. Deep H₂ flux estimated by Prinzhofer, assuming the concentration in the deep venting gas is 50% H₂ compared to air flux in and out of vent that is driven by tidal atmospheric pressure changes. The air flux exceeds the venting flux for venting fluxes at the low end of the range estimated by Prinzhofer for a subsurface permeability of 1 Darcy.

Estimation Method	H2 or Air Flux
Prinzhofer Vdeep for 50% H2 concentration	0.05 to 1.2 m3 m−2 d−1
Calculated p-driven influx for φ = 0.1 to 0.4	0.079 to 0.159 m3 m−2 d−1
Calculated T-driven influx for φ = 0.1 to 0.4	0.005 to 0.018 m3 m−2 d−1

.

The calculations assume a subsurface permeability of 1 Darcy, and porosities between 0.1 and 0.4. The two harmonic components of pressure-driven gas flow in Table 1 are summed to produce the values in the Table 2. The pressure-driven gas flux is not very sensitive to porosity (factor of 2 change for a factor of 4 change in porosity). This $\sqrt{\varphi }$ dependence arises because ${\delta }_{p’}^{}$ depends on $1/\sqrt{\phi }$ (see Table 1). For a subsurface permeability of 10 darcies, ${\delta }_{p’}^{1d}$ would be increased by $\sqrt{10}=3.2$, and the calculated pressure-driven air flux into the ground would range between 0.25 and 0.5 m³ m⁻² d⁻¹.

Diurnal temperature changes penetrate so little into the subsurface that they cannot affect H₂ a meter below the surface. From Table 1, the thermal skin depth for daily temperature variations is at most 10 cm. This means that at 20 cm depth the temperature variation will be 10% of that at the surface, and at 40 cm, 1%. Temperature variations thus do not appear to be a viable way to explain [H₂] variations measured at 1 m depth in the Sao Francisco basin. We therefore do not consider them further in this paper.

2.3.2. Pressure-driven Subsurface Gas Fluxes

Figure 3 shows the changes in pressure, the negative of the rate of pressure change (a proxy for decompression since gas expansion occurs when the pressure decreases), and the gas venting rate, all as a function of depth in meters into the subsurface. The calculated pressure changes take into account the resistance to gas venting offered by a 1-Darcy subsurface with 20% porosity. Note that at depth the pressure changes are out of phase with the surface pressure changes. The venting flux integrates the changes in gas volume from a depth where they are negligible to the surface. The profiles are integrated here from $\overline{z}$ = −10 to 0 (z = −394 m to 0) using Equation (14) with V_deep = 0, but the profiles are displayed in Figure 3 only to 120 m depth. The air flux profiles in the last panel are labeled “cyclic” to emphasize that they are produced by the cyclic daily changes in atmospheric tidal pressure only.

Figure 4A compares the calculated surface venting rate (black curve), the negative of the surface rate of pressure change (blue curve), and the H₂ concentration measured by Sensor 6 for 14 August 2018 (orange curve). The atmospheric pressure is decreasing at a maximum rate at 13:00 (blue curve peak). The maximum H₂ concentration measured by Sensor 6 (and most of the other sensors show in Figure 1) occurs at the same time (orange curve peak). On the other hand, the surface venting rate (black curve) peaks two hours later at 15:00.

Figure 4B shows the relationship between atmospheric pressure changes and the rate of surface gas efflux (V_cyclic). There is considerable hysteresis between the gas flux at 1.2 m depth and the surface pressure changes that produce this flux. The changes in gas flux values are not the same when the surface pressure is increasing as they are when it is decreasing. Figure 4C shows the relationship between the measured [H₂] at Sensor 6 and the calculated cyclic gas flux at 1 m depth. Again, there is substantial hysteresis, and the hysteresis is very similar to that observed and shown in Figure 2C.

Figure 5 shows that if the cyclic (compression/decompression) gas flux comes from only the shallowest parts of the full system (the first 10 m of the curves shown in Figure 3), the phase shift with respect to the peak [H₂]_Sensor-6 is reduced. If the compression/decompression gas flux derives from ~10 m (= 0.25 ${\delta }_{p’}^{1d}$, where ${\delta }_{p’}^{1d}$ = 39.4 m)of the subsurface (${\overline{z}}_b$ = 0.25 in Equation (9)), the maximum in the surface venting rate (gray curve) nearly coincides in time with the maximum measured H₂ concentration. Note, however, that the magnitude of the gas flux is strongly reduced from the black curve (Figure 5A) where gas is expelled from all depths.

Figure 5B shows that the air flux profiles are nearly linear if they originate in only the upper 10 m of the profiles shown in Figure 3. Their linearity indicates that the upper part of the vent is behaving like the simple box model defined in Equation (16) where it is assumed that all the air in the subsurface is compressed equally and at the same time as the surface pressure changes. In the box case, the integrated flux will increase linearly from zero at ${\overline{z}}_b$, as is nearly the case in Figure 5B.

2.3.3. Box Model Advection/Diffusion Calculations of [H₂] at 1 m Depth

For a rise in atmospheric pressure to cause a drop in H₂ concentration at 1 m depth, air must penetrate at least that deep. This is a significant constraint: The air velocity v_air equals V_air/φ, so for φ = 0.2 and V_air = 0.05 m d⁻¹ (the low flux estimate of Prinzhofer), the air velocity is 0.25 m d⁻¹. The atmospheric tide drives air into the ground for at most half a day, so the depth of air penetration for this low H₂ flux estimate is at most ~12 cm, which is much less than the depth of the sensors. To reach 1 m depth, the influx must be at least 8 times greater, e.g., V_air~0.4 m d⁻¹. Similarly, dilution of the effluent H₂ gases by influent air must be quite substantial for the measured [H₂]_1m concentrations to be 100 to 1000 ppm. A simple stirred tank mixing model suggests reducing [H₂]_1m by a factor of 1000 (from 50% H₂ to 500 ppm) would require mixing one volume of deep hydrogen-rich gas with ~7 volumes of air. Thus, from both these perspectives, the air influx must be about an order of magnitude greater than Prinzhofer’s low estimate or ~0.05 m d⁻¹.

The advection–diffusion Equation (17) can be used to model the [H₂] concentration at 1 m depth. The thickness of the “box” must be ~1000 m to produce cyclic air velocities large enough to dilute the hydrogen concentration at 1 m depth by the large amounts observed. This is an unrealistically large depth but suitable for our heuristic calculations here. The calculations start with the initial steady state [H₂] profile in b_box for the assumed deep gas flux V_deep. This initial profile is then modified by the cyclic tidal air movements. The gas flux is constant but time varying at its 1-m-deep box value over the depth interval of analysis. The concentration changes produced by the cyclic variations in air flow are calculated to a much shallower depth, b<<b_box, using Galerkin finite element methods (Baker and Pepper, 1991 [26]) because the changes in [H₂] induced by the cyclic air flux do not extend very far into the subsurface (<25 m).

The calculations proceed in 1 h timesteps with 100 sub-timesteps. For plotting, [H₂] is normalized by dividing by the deep H₂ concentration: $[{\overline{H}}_2]=[H_2]/{[H_2]}_{deep}$. The bottom boundary condition is Neuman wherein the H₂ flux equals $V_{deep}{[{\overline{H}}_2]}_{deep}$ and ${[{\overline{H}}_2]}_{deep}$ = 1 (100% H₂). The top boundary is Dirichlet, ${[{\overline{H}}_2]}_{z=0}$ = 0. Four daily pressure cycles are computed, and the last selected for plotting. b_box = 1000 m and V_deep = 0.1 m d⁻¹. There are 100 finite elements in the 25 m interval calculated. The surface air venting rate varies from −3.4 to + 5 m d⁻¹. Calculating for 200 instead of 100 finite elements in the interval from the surface to b = 25 m makes no difference to the results. The 200 element curves exactly overlie the 100 element blue curves, showing the 100 element calculation has more than enough depth resolution.

Figure 6 shows the changes in the $[{\overline{H}}_2]$ calculated as a function of depth and time of day for the atmospheric pressure variation of 14 August 2018 in the Sao Francisco Basin of Brazil. The insert in each panel shows the gas flux at the surface over the interval plotted. It can be seen that the tidal air flux produces dramatic changes in near surface hydrogen concentrations. Periods of air inflow drive the near-surface ${\left[{\overline{H}}_2\right]}_{1m}$ concentration to lower values. Periods of air outflow increase the near-surface ${\left[{\overline{H}}_2\right]}_{1m}$. The H₂ concentration varies between the surface and ~20 m depth. Below this depth the H₂ concentration is at the deep, up-flow value.

Figure 7 shows the calculated hydrogen concentration at 1 m depth, ${\left[{\overline{H}}_2\right]}_{1m}$ on 14 August 2018. Periods of air inflow drive the near-surface H₂ concentration to very low values except for a small peak at ~3 AM. Periods of venting elevate the 1 m depth concentrations to 37% of the deep input concentration. The calculated variations in ${[H_2]}_{1m}$ are offset from those observed (orange curve) by about ~3 h. The calculated modulation of ${\left[{\overline{H}}_2\right]}_{1m}$ is strongly dependent on the vertical extent of the high permeability “box”. The dashed blue line shows the comparatively small variation in ${\left[{\overline{H}}_2\right]}_{1m}$ calculated for b_box = 200 m.

Finally, Figure 8 shows that, as expected, there is no substantial hysteresis between ${\left[H_2\right]}_{1m}$ and atmospheric pressure for the box model calculations. This is in strong contrast to observations (Figure 2C). Also, the calculated ${\left[H_2\right]}_{1m}$ concentration continues to increase until the gas efflux reverses (e.g., it is maximum at 17:00), rather than peaking when the rate of pressure change is greatest at 13:00 (red circle).

3. Discussion

From the above analysis, it is clear that the atmospheric pressure tide could cause the pulsation in hydrogen concentration measured at 1 m depth. The main evidence supporting the hypothesis that atmospheric pressure changes are modulating the measured hydrogen concentration is the hysteresis in the observed Δp’ vs. [H₂]_1m curve shown in Figure 2C. Figure 4C shows that the slow calculated diffusion of pressure into the subsurface produces a hysteresis between the rate of venting at 1.2 m depth, V_1.2m, and the measured [H₂]_1m that is very similar in form to that observed between Δp’ and [H₂]_1m. (Note, the circulation is the same if Δp’ is replaced by −Δp’.) If [H₂]_1m is proportional to V_1.2m, this hysteresis similarity strongly suggests the diffusion of pressure into the subsurface is the cause of the measured pulsing of the hydrogen venting.

A substantial reservoir of gas (compared to the volume of gas in the vents) must be compressed or expanded by the atmospheric pressure changes for hydrogen-free atmospheric air to be drawn to sensor depth. For the simple “box” model calculated above, the box must be ~1000 m thick to change [H₂]_1m in a fashion similar to that observed (Figure 7). Instantaneous pressure transmission to 1000 m depth would require an unrealistically high subsurface permeability, so the box depth of 1000 m simply indicates that the reservoir gas volume affected by atmospheric pressure changes must be at least ~1000 times larger than the gas volume between the surface and the H₂ sensor at 1 m depth. Pressure wave calculations show that, in addition, the volume of gas accessible to the pressure wave must be about 25% of the full volume with which it could interact. This is required for the maximum venting rate, V, to coincide with the maximum [H₂]_1m (Figure 5A). If the vent has a very low gas volume compared to the reservoir with which it interacts, there will be very little transit delay for incoming air to reach the H₂ gas sensors at 1 m depth. It is important to emphasize that the box modeling is 1D. Flow arises from vertical compression and decompression only. In reality gas would be supplied to vents laterally as well as vertically. Thus, the 25% of the potential draw should be interpreted as 25% of the 3D volume that feeds a particular vent.

It is reasonable that [H₂]_1m should be maximum at the maximum venting rate at 1 m depth. The advection–diffusion solution shown in Figure 8 seems to suggest differently. It shows the maximum [H₂]_1m occurs at the end of venting just before inflow begins. However, it is a 1D calculation that considers only vertical diffusion. In actuality, H₂ diffuses laterally from the vent, and, as gas rich in H₂ approaches the surface H₂ will diffuse laterally and be diluted. This dilution will be minimum when the gas efflux is maximum, and thus the maximum [H₂]_1m should coincide with the maximum venting rate.

Observations as well as these modeling results suggest that hydrogen is venting from a shallow reservoir lying between the surface and the water table under the barren zone, and that the venting occurs mainly on the periphery, as shown in Figure 9. The barren zone is a small topographic depression that fills periodically with water. It is plausible that the top of the barren zone could be less permeable than its periphery because, due to periodic flooding, it receives more fine sediment deposition, has more evaporative salt deposition, and is more altered. Because slumping permeability might also be concentrated at the barren zone margins. If the upper layer of the central portion of the barren zone is relatively impermeable, but underlain by permeable sediments, a sealed H₂ reservoir could exist in the permeable sediments between the surface and the water table. The hypothetical reservoir could extend outside the barren zone if the sediments above the water table were as permeable outside as inside the barren zone. The reservoir in Figure 9 is shown being filled from depth by relatively pure (50 to 100%) hydrogen gas. Gas pushed into and out of the reservoir by atmospheric pressure tides through the periphery vents dilutes the reservoir near these vents as illustrated by the green hydrogen concentration contours surrounding the leftmost vent in Figure 9. Different perimeter vents will interact with the reservoir in slightly (and perhaps substantially) different ways if the permeability and porosity vary around the periphery of the reservoir. The time and concentration of the peak hydrogen concentration at different sensors could therefore differ as observed. The vents will operate as observed provided the three-to-five-meter-thick reservoir constitutes ~25% of the potential pressure wave penetration depth and the gas volume in the vents is a very small fraction of the volume compressed and decompressed by the atmospheric pressure tides impacting each vent.

The possibility that the vent system is operating as illustrated in Figure 9 can be tested in several ways: A gas probe in the center of the barren zone would test if there is a gas reservoir between a sealed surface and the water table. The hydrogen concentration in the center of the barren zone should be >50% (or at least much greater than near the vents). Gas probes into the reservoir near sites of venting on the periphery could show how [H₂] varies away from the vents. The gradient in [H₂] and pressure variations at these probes could confirm the hypothesis that atmospheric-pressure variations cause the observed changes in measured H₂. Measurements of permeability and porosity would also test this hypothesis and would provide data for the kind of 3D finite element analysis that will be needed to accurately model the H₂ venting. Drill holes outside the barren zone would test the extent of the H₂ reservoir.

There is much that is not covered in our analysis. For example, as the water table at the base of our hypothetical H₂ reservoir rises and falls, accumulated H₂ will be expelled and diluted. Tracking these changes could be important to the H₂ content of the reservoir. The magnitude of H₂ venting is best provided, at least in the short term, by integrating the H₂ efflux from the periphery of the barren zone as has been done by Prinzhofer et al. (2019). We add nothing to Prinzhofer’s estimates of the total H₂ venting rate in this paper. Rather, the analysis presented in this paper suggests the kind of system that could operate as observed. Ultimately 3D finite element modeling will be needed to define the hydrogen resource. For 3D modeling to contribute beyond the analysis offered here, however, more needs to be known about the shallow H₂ reservoir and its relation to the vents on the periphery of the barren zone. The needed information can be obtained by gas probe or shallow drilling.

4. Conclusions

This paper analyzes whether pressure- or temperature-driven air flow can explain the temporal variations in hydrogen concentration measured at 1 m depth along the perimeter of a 550 m diameter, largely barren depression in the Sao Francisco Basin in Brazil. Although the temporal variations could be caused by other processes such as solid earth tides or temperature-dependent bacterial H₂ generation, etc., we find:

• The variations in hydrogen concentration measured at 1 m depth could be caused by propagation of a pressure wave into the subsurface, but not by the propagation of a thermal wave. Diurnal temperature changes penetrate less than a meter into the subsurface and produce only weak perturbations of gas flow above the one-meter sensor depth.
• The propagation of a pressure wave truncated at 25% of its potential penetration could produce changes in hydrogen concentration at 1 m with a phase shift relative to atmospheric tidal pressure changes similar to that most commonly observed.
• To change [H₂] concentrations at 1 m depth, the gas volume in each vent needs to be <1/1000th of the reservoir gas volume with which atmospheric pressure variations interact.
• The venting system we infer here from observations and calculations is illustrated in Figure 9. The characteristics of this hypothetical system should be tested by measuring reservoir hydrogen concentrations with gas probes or by drilling as indicated in the discussion section.