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

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 H2 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 H2 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.


Introduction
Hydrogen gas may become an important part of the zero-carbon energy economy because its combustion produces only water vapor and no CO 2 . 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.
The most important parameter for gauging the magnitude of the hydrogen resource is the magnitude of the hydrogen (H 2 ) 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 2 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 2 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 2 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 2 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 2 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 3 H 2 per day, based on the venting at Sensor 9 ( Figure 1) where the volume of H 2 venting is greatest and the concentration of H 2 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 4 + decomposition at depths where the temperature exceeds 200 to 250 • C (Guélard et al., 2018) [22]. The H 2 concentration was found to be 50% to 80% (with the rest N 2 and several percent CH 4 ) in a deep drill hole in the Sao Francisco basin, and measurements at other sites in the basin found a source H 2 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 2 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 2 -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 2 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.
Abarren. The study area is the large barren area shown in Figure 1 with Abarren = 300,000 m 2 . It is venting hydrogen at an estimated rate of between 7000 and 178,000 m 3 H2 per day (Prinzhofer et al., 2019) [11]. Assuming the concentration of the source is 50% H2, the venting rate at depth, QH2, is between 14,000 and 356,000 m 3 per day. If the flow is uniform across the vent area, the gas flux is Vdeep = QH2/Abarren = 0.05 to 1.2 m d −1 . The average concentration of H2 in the discharge measured at ~1 m depth is ~100 ppm (Prinzhofer et

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.

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 [H2] concentrations measured at 1 m depth by Sensor 6 on the same day. Peak [H2]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).

The Study Site
As illustrated in Figure 1, H 2 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 2 . It is venting hydrogen at an estimated rate of between 7000 and 178,000 m 3 H 2 per day (Prinzhofer et al., 2019) [11]. Assuming the concentration of the source is 50% H 2 , the venting rate at depth, Q H2 , is between 14,000 and 356,000 m 3 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 −1 . The average concentration of H 2 in the discharge measured at~1 m depth is~100 ppm (Prinzhofer et

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.

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 2 ] concentrations measured at 1 m depth by Sensor 6 on the same day. Peak [H 2 ] 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). 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 [H2] is measured at 13:00 when the surface pressure is decreasing at its maximum rate ( 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 Geosciences 2020, 10, 149 5 of 18 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 2 ] 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.

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: Conservation of mass requires: 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 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: where Because we want the rate of venting, we take the time derivative of (5): Defining z = z/δ α , we can integrate (5) from some depth 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: 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 z simply indicates that it is an integration variable. If the depth of integration, z b , is sufficiently large that the underlying pressure or temperature variations are negligible, for example z b = −10, the normalized integral, 0 z=−10 e z sin(2πt + 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 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 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: Since the gas flux is linear, we can add the deep H 2 -bearing gas flux, V deep , to (6) to obtain the total gas flux, V total Z=0 (t): The H 2 emissions could be pulsating if: Ground changes > atmospheric.
β T -Thermal expansion of air       [11], assuming uniform venting in barren zone.
If condition (12) is satisfied, the T-driven or p'-driven flow into the ground will exceed the base leakage rate of deep H 2 -bearing gas (V deep ), H 2 -free air could periodically surround the sensors, and the sensors could then detect low H 2 concentrations. Since V deep = 0.05 to 1.2 m 3 m −2 d −1 , 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: Geosciences 2020, 10, 149 9 of 18 For ϕ = 0.2, b box = 1000 m, β p = 10 −5 Pa −1 , and − ∂p ∂t = 25 mb d −1 = 2500 Pa d −1 , the gas flux from the box V box = 5 m d −1 . 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 2 (expressed as the volume fraction of H 2 in air) by diffusion and advection is: where [H 2 ] is the volume fraction of H 2 in air. By mass conservation: where D E is the effective diffusion constant of H 2 in the vent, and V air is the vertical flux of air (m 3 of air passing through a plan area or 1 m 2 per second) in the portion of the vent considered (in this case the uppermost part). The effective diffusion constant D E = Dφ τ , where D is the diffusion constant of H 2 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 2 ] non-dimensional by defining: Equation (15) becomes: where N pe is the Peclet number (ratio of advection to diffusion).

The Minimum Criterion for [H 2 ] Station-6 Modulation
Whether the pulsing H 2 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 2 ] venting if the venting were at the low end of Prinzhofer's estimated range. For convenience we show the comparison in Table 2. Table 2. Deep H 2 flux estimated by Prinzhofer, assuming the concentration in the deep venting gas is 50% H 2 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 H 2 or Air Flux
Prinzhofer V deep for 50% H 2 concentration 0.05 to 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 φ dependence arises because δ p depends on 1/ √ ϕ (see Table 1). For a subsurface permeability of 10 darcies, δ 1d p would be increased by Diurnal temperature changes penetrate so little into the subsurface that they cannot affect H 2 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 2 ] variations measured at 1 m depth in the Sao Francisco basin. We therefore do not consider them further in this paper. 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 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 H2 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 H2 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  Figure 4A compares the calculated surface venting rate (black curve), the negative of the surface rate of pressure change (blue curve), and the H 2 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 2 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.

Pressure-driven Subsurface Gas Fluxes
shown in the legend in the last panel. The curves are for different hours of the day as indicated by the color key at the upper right of each diagram. Figure 4A compares the calculated surface venting rate (black curve), the negative of the surface rate of pressure change (blue curve), and the H2 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 H2 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 (Vcyclic). 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 [H2] 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 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 2 ] 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 2 ] Sensor-6 is reduced. If the compression/decompression gas flux derives from~10 m (= 0.25 δ 1d p , where δ 1d p = 39.4 m)of the subsurface (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 2 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 z b , as is nearly the case in Figure 5B. 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 b z , as is nearly the case in Figure 5B.

Box Model Advection/Diffusion Calculations of [H2] at 1 m Depth
For a rise in atmospheric pressure to cause a drop in H2 concentration at 1 m depth, air must penetrate at least that deep. This is a significant constraint: The air velocity vair equals Vair/φ, so for φ = 0.2 and Vair = 0.05 m d −1 (the low flux estimate of Prinzhofer), the air velocity is 0.25 m d −1 . The atmospheric tide drives air into the ground for at most half a day, so the depth of air penetration for this low H2 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., Vair~0.4 m d −1 . Similarly, dilution of the effluent H2 gases by influent air must be quite substantial for the measured [H2]1m concentrations to be 100 to 1000 ppm. A simple stirred tank mixing model suggests reducing [H2]1m by a factor of 1000 (from 50% H2 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 −1 .
The advection-diffusion Equation (17) can be used to model the [H2] 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 [H2] profile in bbox for the assumed deep gas flux Vdeep. 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<<bbox, using Galerkin finite element methods (Baker and

Box Model Advection/Diffusion Calculations of [H 2 ] at 1 m Depth
For a rise in atmospheric pressure to cause a drop in H 2 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 −1 (the low flux estimate of Prinzhofer), the air velocity is 0.25 m d −1 . The atmospheric tide drives air into the ground for at most half a day, so the depth of air penetration for this low H 2 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 −1 . Similarly, dilution of the effluent H 2 gases by influent air must be quite substantial for the measured [H 2 ] 1m concentrations to be 100 to 1000 ppm. A simple stirred tank mixing model suggests reducing [H 2 ] 1m by a factor of 1000 (from 50% H 2 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 −1 .
The advection-diffusion Equation (17) can be used to model the [H 2 ] 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 2 ] 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 [ into the subsurface (< 25 m).
The calculations proceed in 1 h timesteps with 100 sub-timesteps.        Figure 8 shows that, as expected, there is no substantial hysteresis between [H 2 ] 1m and atmospheric pressure for the box model calculations. This is in strong contrast to observations ( Figure 2C). Also, the calculated [H 2 ] 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). it is maximum at 17:00), rather than peaking when the rate of pressure change is greatest at 13:00 (red circle).   it is maximum at 17:00), rather than peaking when the rate of pressure change is greatest at 13:00 (red circle).

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 2 ] 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 2 ] 1m that is very similar in form to that observed between ∆p' and [H 2 ] 1m . (Note, the circulation is the same if ∆p' is replaced by −∆p'.) If [H 2 ] 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 2 ] 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 2 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 2 ] 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 2 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 2 ] 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 2 ] 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 2 diffuses laterally from the vent, and, as gas rich in H 2 approaches the surface H 2 will diffuse laterally and be diluted. This dilution will be minimum when the gas efflux is maximum, and thus the maximum [H 2 ] 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 2 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.
advection-diffusion solution shown in Figure 8 seems to suggest differently. It shows the maximum [H2]1m occurs at the end of venting just before inflow begins. However, it is a 1D calculation that considers only vertical diffusion. In actuality, H2 diffuses laterally from the vent, and, as gas rich in H2 approaches the surface H2 will diffuse laterally and be diluted. This dilution will be minimum when the gas efflux is maximum, and thus the maximum [H2]1m should coincide with the maximum venting rate. Figure 9. Schematic of H2 vent system suggested by the modeling results. Figure 9. Schematic of H2 vent system suggested by the modeling results.
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 2 ] varies away from the vents. The gradient in [H 2 ] and pressure variations at these probes could confirm the hypothesis that atmospheric-pressure variations cause the observed changes in measured H 2 . 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 2 venting. Drill holes outside the barren zone would test the extent of the H 2 reservoir.
There is much that is not covered in our analysis. For example, as the water table at the base of our hypothetical H 2 reservoir rises and falls, accumulated H 2 will be expelled and diluted. Tracking these changes could be important to the H 2 content of the reservoir. The magnitude of H 2 venting is best provided, at least in the short term, by integrating the H 2 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 2 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 2 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.

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 2 generation, etc., we find:

1.
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.

2.
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.

3.
To change [H 2 ] 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. 4.
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.