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 *z* = −10 to 0 (z = −394 m to 0) using Equation (14) with *Vdeep* = 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 3.** Calculated profiles of (**left**) pressure, (**middle**) the negative of the rate of pressure change (a proxy for gas decompression), and (**right**) the gas venting rate. Changes are caused by the changes atmospheric pressure shown in Figure 2. The subsurface permeability is 1 Darcy; the porosity 0.2 as 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 *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 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 4.** (**A**) The calculated gas venting rate at the surface (black curve) plotted versus local time for 14 August 2018. The blue curve shows the negative of the rate of change in atmospheric pressure at the surface. The orange curve shows the H2 concentration measurements in Sensor 6 as a function of local time. The peaks in the H2 concentration and –dp/dt coincide at 13:00, but the peak in surface venting rate is shifted by 2 h and occurs at 15:00. (**B**) The calculated gas flux Vcyclic at the surface is phase lagged relative to atmospheric pressure variations. (**C**) Calculated Vcyclic at 1.2 m depth is phase lagged relative to measured [H2]1m at Sensor 6 is a fashion similar to atmospheric pressure variations shown in Figure 2C. The red circles indicate the time maximum [H2]1m is observed in Sensor 6.

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 [*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 δ1*dp* , where δ1*dp* = 39.4 m)of the subsurface (*zb* = 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 *zb*, as is nearly the case in Figure 5B.

**Figure 5.** (**A**) The maximum cyclic gas flux calculated from a depth of 0.25δ (10 m depth for δ = 39.4 m) is more in phase with the maximum measured H2. (**B**) Depth profiles of gas flux from 09:00 to 18:00. The depth profiles are almost linear and therefore approximate the flux profiles that would occur in a simple box model of the kind described in the box model section.

### 2.3.3. Box Model Advection/Diffusion Calculations of [H2] 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 *vair* equals *Vair*/ϕ, so for ϕ = 0.2 and *Vair* = 0.05 m d−<sup>1</sup> (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., *Vair*~0.4 m d−1. Similarly, dilution of the effluent *H*2 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% *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 *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 Pepper, 1991 [26]) because the changes in [*H*2] 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*2] is normalized by dividing by the deep *H*2 concentration: [*<sup>H</sup>*2]=[*<sup>H</sup>*2]/[*<sup>H</sup>*2]*deep*. The bottom boundary condition is Neuman wherein the *H*2 flux equals *Vdeep*[*<sup>H</sup>*2]*deep* and [*<sup>H</sup>*2]*deep* = 1 (100% H2). The top boundary is Dirichlet, [*<sup>H</sup>*2]*<sup>z</sup>*=<sup>0</sup> = 0. Four daily pressure cycles are computed, and the last selected for plotting. *bbox* = 1000 m and *Vdeep* = 0.1 m d−1. There are 100 finite elements in the 25 m interval calculated. The surface air venting rate varies from −3.4 to + 5md−1. 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 [*<sup>H</sup>*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 *<sup>H</sup>*2<sup>1</sup>*m* concentration to lower values. Periods of air outflow increase the near-surface *<sup>H</sup>*2<sup>1</sup>*m*. The *H*2 concentration varies between the surface and ~20 m depth. Below this depth the *H*2 concentration is at the deep, up-flow value.

**Figure 6.** Calculated profiles of *<sup>H</sup>*2for bbox = 1000 m and Vdeep = 0.1 m d−<sup>1</sup> for various times during 14 August 2018. The insert in each panel shows the calculated venting rate, highlighting the interval of the day when the profiles are plotted. The legend lists the time of day of each profile.

Figure 7 shows the calculated hydrogen concentration at 1 m depth, *<sup>H</sup>*2<sup>1</sup>*m* on 14 August 2018. Periods of air inflow drive the near-surface *H*2 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 [*<sup>H</sup>*2]<sup>1</sup>*m* are offset from those observed (orange curve) by about ~3 h. The calculated modulation of *<sup>H</sup>*2<sup>1</sup>*m* is strongly dependent on the vertical extent of the high permeability "box". The dashed blue line shows the comparatively small variation in *<sup>H</sup>*2<sup>1</sup>*m* calculated for *bbox* = 200 m.

Finally, Figure 8 shows that, as expected, there is no substantial hysteresis between [*<sup>H</sup>*2]<sup>1</sup>*m* and atmospheric pressure for the box model calculations. This is in strong contrast to observations (Figure 2C). Also, the calculated [*<sup>H</sup>*2]<sup>1</sup>*m* 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).

**Figure 7.** Computed *<sup>H</sup>*2<sup>1</sup>*m* for *bbox* = 1000 m and *Vdeep* = 0.1 m d−<sup>1</sup> (blue curve) is compared to the ppm [*<sup>H</sup>*2]<sup>1</sup>*m* observed (orange curve). The dashed blue curve the H2 profile computed for *bbox* = 200 m.

**Figure 8.** Phase relationship between [*<sup>H</sup>*2]<sup>1</sup>*m* and the atmospheric pressure changes that drive the cyclic gas flow calculated in the box model. Unlike observations (Figure 2C), no substantial hysteresis is suggested by the box model calculations. The red circle marks the time maximum [*<sup>H</sup>*2]<sup>1</sup>*m* is observed in Sensor 6.
