© A. Myagkiy et al., Published by EDP Sciences 2020

1 Introduction

1.1 H₂ transfer in subsurface

Subsurface H₂ gas concentrations far above atmospheric level and subsequent H₂ flow towards the surface have been reported in an increasing number of publications. H₂ gas concentrations at around 1 m below the surface have been measured in Russia and USA among others (Larin et al., 2014; Zgonnik et al., 2015) and even monitored in Brazil (Prinzhofer et al., 2019) on a time-resolved basis. The source of this hydrogen is still debated since different H₂-forming reactions and processes are likely to produce H₂ in the subsurface (Guélard et al., 2017). The hydrothermal oxidation by the seawater of ultramafic rocks exhumed along the mid oceanic ridges is one of the best-known redox process that produces H₂ under abiotic conditions (e.g., Brunet, 2019 and references therein). H₂ fluxes related to this phenomenon, called serpentinization, have been observed along the mid-Atlantic ridge Charlou et al. (2002, 2010) or the Pacific one (Welhan and Craig, 1979) where it takes place at high temperatures under hydrothermal conditions. H₂ emissions are also observed in convergence area where the oceanic lithosphere outcrops onshore such as in Oman, Philippines or New Caledonia (Neal and Stanger, 1983; Abrajano et al., 1988; Sano et al., 1993; Monnin et al., 2014; Deville and Prinzhofer, 2016; Vacquand et al., 2018). The same type of process based on the interaction between ultramafic rocks and water may occur at lower temperature and the fluid could be meteoric water (e.g., Mayhew et al., 2013; Miller et al., 2017).

Similar reactions of oxidation/hydration may also take place in continental settings if, for instance, Fe²⁺ rich rocks are present, as it has been suggested as a source for the continental H₂ accumulations in Mali (Prinzhofer et al., 2018). Radiolytic H₂ production in the Precambrian continental crust, which accounts for 72% of the continents, is another possibility (Sherwood et al., 2014) with a contribution which, for these last authors, could be comparable to the one of the oceanic domains, around 10¹¹ moles per year. In addition, H₂ may originate from mantle degassing (Larin et al., 2010) since H₂ is the initial and major component of the solar system. Mechanoradical H₂ generation along fault planes when they are creeping, or active, have also been observed (Klein et al., 2020); following these authors it can occur wherever there is comminution of silicate-bearing rocks. A frictional process has been proposed to explain the change in H₂/CH₄ ratio in a gas field in the Jura front belt (Deronzier and Giouse, 2020) however its overall role for large scale H₂ production remains unknown. Bacterial activity is also prone to produce H₂ in specific environments (Conrad and Seiler, 1980), however this alternative has also not yet been quantified at the global scale.

Without attempting to discriminate between these possible sources, the present work focuses on the H₂ transfer in continental settings characterized by a soil cover and which leads to the formation of fairy circles. In the absence of soil cover, H₂ emanations are usually reported along faults. In the Oman ophiolite, for example, H₂-rich gases have the same exit points as alkaline spring waters, which are located along fault and shear zones (Neal and Stanger, 1983). In such a case, water and gases react immediately with the atmosphere resulting in carbonate deposits such as in Oman (Noel et al., 2018) and in New Caledonia (Deville and Prinzhofer, 2016). In case of water accumulation (river or lake), H₂ being poorly soluble in water at low pressure (Lopez-Lazaro et al., 2019 and reference therein), continuous degassing can be observed in the form of bubbles. They show a continuous but not permanent flow of H₂.

After a synthesis of the characteristics of H₂ emitting structures in terms of shape and H₂ concentrations in the soils, we explore, in the current study, through two-dimensional reaction-transport modeling, the situation where the H₂ migration toward the surface ends in a soil cover. Our model of H₂ subsurface distribution aims at accounting for the H₂ concentrations described as non-constant neither in time nor space reported worldwide in the circular-like shape structures.

1.2 Shape of the H₂ emitting structures in the presence of a soil layer

Due to their circular shape, the structures from which H₂ is leaking out to the atmosphere are often called “fairy circle”. We will use this wording even if the structures are not perfectly circular, and even sometimes perfectly elliptic as in South Carolina (Zgonnik et al., 2015). Figure 1 shows some examples of these structures.

Larin et al. (2014) have investigated a large number of them in the Voronezh Oblast, about 600 km south-east of Moscow. The size distribution deduced from the 562 measured circles shows variability from 3 km down to 100 m (Fig. 2 in Larin et al., 2014). Small structures are more frequent, and the authors noticed that above 1.5 km there is no more relationship between size and frequency. The authors noted that if the center of the structures corresponds sometimes to a true depression, up to 5 m deep, often there is no real depression but just changes in the vegetation. When the center corresponds to a depression, wetlands are present. Small lakes and changes in the soil chemistry are noted, Larin et al. (2014) noticed that soils are sandier at the circle rim and more shaly in the center. Similarly, in South Carolina (USA), hundreds of depressions have been also reported, first by researchers working on wetlands and biodiversity (Semlitsch, 2000) and then by those studying H₂ emissions (Zgonnik et al., 2015). They noted, as in Russia, a large spectrum of sizes from 8 km down to 100 m for the major axis of the depressions which, in that case, are ellipses (Fig. 1D).

Larin et al. (2014) reported a structural alignment and suggested that the fairy circles are located above faults from where the H₂ leaks out. The crucial role of faults for fluid (liquid as gas) migration is well established and will not be discussed here. The general view is that pervasive circulation of the fluids in the geological medium is large when bulk porosity and permeability are large – that is the case in all carrier beds and aquifers – whereas as soon as the permeability is low (shale, metamorphic rock, granite …) the fluids circulate mainly in the fractures and faults. The cross-section (Fig. 3 in Larin et al., 2014) shows that, in their case study, the Precambrian basement is located between 0 and 500 m depth. The fact that H₂ escapes from and along fractures from the basement is thus a possibility, but when reaching the overlying soil cover, upward transport will proceed within the soil medium which has usually a relatively “good” porosity and permeability that are at least few orders of magnitude higher than that of the Precambrian basement.

When there is no soil cover, as in Oman, New Caledonia (Deville and Prinzhofer, 2016; Vacquand et al., 2018) but also Turkey in Chimeres area or in Philippines, the gas mainly seeps from fractures, these settings will not be discussed here. Roughly speaking the shape of these seeps is the shape of the fracture network and is also out of the scope of this study.

Fig. 1 Synthesis of the fairy circle emissions characteristics. A. A fairy circle in Brazil showing how the vegetation is affected. B. Variability on time over 5 weeks for a given sensor in the structure monitored in Brazil. See Prinzhofer et al. (2019) and Myagkiy et al. (2020) for details about that case. C. Daily variations for 6 sensors installed on the structure monitored in Brazil showing the 24 h period. D. Variability in size, of the emitting structures in USA (modified from Prinzhofer and Deville, 2015). E. Two profiles across a structure showing the spatial variability of H2 concentration value between the center and the rim of the structure, Lake Podowoye, Russia. These values have not been taken at the same time but over 2 days. Modified from Larin et al. (2014). F. Profile showing an example of [H2] concentration variation versus depth in USA (modified from Zgonnik et al., 2015).

1.3 Time-resolved H₂-concentration measurements within the fairy circles

Up to now, even though the list of basins where H₂ emanations have been reported is getting longer (Zgonnik, 2020), both time-resolved and spatial distribution measurements of H₂ concentration in soils from fairy circles are still rare. However, the datasets cited in the previous section all show that the soil H₂-concentration is variable across the structure and that the maximum concentration is barely, if not never, in the center of the structure (Figs. 5 and 6, Larin et al., 2014; Prinzhofer et al., 2019).

The only published permanent monitoring of an H₂ concentration in the soil porosity of an emitting structure has been performed in the São Francisco basin in Brazil (Prinzhofer et al., 2019). A few other examples of published permanent H₂ monitoring data are either related to seismic risk and were carried out in the vicinity of active faults, as in San Andreas (Sato et al., 1986) and in Japan or for safety control during mining operations (Kola Peninsula, Russia, Firstov and Shirokov, 2005; Nivin et al., 2016).

In Brazil, the monitored circular structure is a depression with a diameter close to 500 m. The published H₂-concentration dataset covers more or less one and half month of record. It shows a daily wavelength consistent with a pulsed source of a still unclear origin. Almost all the sensors which are near the rim of the structure and installed at 80 cm depth, show a maximum at around noon and are recording H₂ flow for about 6 h a day. H₂ soil concentrations are measured each hour and were around 150 ppm during the week of August 2019 that is reported in Prinzhofer et al. (2019). Another sensor recorded higher values (up to 1200 ppm), it showed also a 24 h period but H₂ concentration maximum was reached during the night. A sensor installed in the center of the structure for few months recorded only a few H₂ concentration pulses but most of the time had a signal close to 0 (Prinzhofer and Engie’s Team personal communication). The 24 h period is not yet fully understood (Myagkiy et al., 2020) nor the gap, or delay, between the highly emitting area and the other ones. However, the fact that the H₂ concentration in the soil, near the surface, is changing with time along the day is shown by all data.

The effect, on H₂ diffusive transport, of bacterial activity has already been investigated on soils from the same structure of the São Francisco basin (Myagkiy et al., 2020). The kinetics of H₂ consumption by the bacteria present in the soil, which was calibrated in laboratory experiments, is proved to be very large and may reduce or even completely hide the deeper H₂ flows. Basically, the consumption of H₂ by bacteria increases with the residence time of H₂ in the soil whatever the transport mode, advective or diffusive.

There are only few published data of H₂ in the soil at various depths. In Oman, Zgonnik et al. (2015) recorded concentrations at 3 different depths of 30, 60 and 110 cm and found concentrations of 150, 375 and 3400 ppm respectively (S54, f1, 2 and 3, Tab. A1 in the cited paper). However, in the same article they report, also data showing an inverse gradient with H₂ concentration which decreases with depth (sample 51 f1 and f2). These data correspond to one-shot measurements without indication of sampling time so even if these variations of concentration with depth are real, they could merely be a snapshot of a situation that may evolve with time.

The monitoring carried out in Brazil has already shown that recording data versus time for a least few days is crucial not only to compute a flow but also even just to know whether H₂ emission do take place. A single zero H₂ concentration (that we will note [H₂] further in this paper) value at a given time cannot be interpreted as a lack of H₂ flow in a dynamic system. Based on a numerical model of H₂ transmission in the soil we will evaluate here which are the key parameters that affect the near surface signal and can account for (1) the size of these structures; (2) the presence of a ring with higher H₂ concentration; (3) the variation through time of the near surface signal. We will therefore discuss the causes of the apparent randomness of the H₂ concentration recorded near the surface.

2 Transport of H₂ across a soil

2.1 Physical governing equations and assumptions

A 2-D numerical advective-diffusive model of gas reactive transport in porous media was developed. For that purpose, the numerical simulator COMSOL Multiphysics v 5.2 was used. All the details concerning the equations and the parameters used in our computations are presented in Appendix A.

Bacterial activity, that was found to have a prominent role in the consumption, or production of gas in the upper part of the subsurface, was incorporated in calculations through the Michaelis–Menten kinetic law of H₂ consumption. The parameters, experimentally derived for the soils of Brazilian H₂-emitting structure, have been used in simulations. Two types of calculation were performed, namely with and without H₂ bacterial consumption, for a comparison. The details of its measurements can be found in Myagkiy et al. (2020). The model assumes a constant gradient in bacteria distribution; the maximum value of bacterial consumption activity is assumed to appear at the surface and reaches zero at the bottom of the model.

It should be noted that the numerical model was resolved in the range of [H₂] at which microbiological uptake kinetics was measured, so up to few thousandths of ppm. For the sake of simplicity, modelled H₂ concentrations in soil were expressed as percentage of the source concentration. The impact of concentration at source in the range of applicability of microbiological consumption law was found to be negligible compared to the other parameters. The [H₂] at the surface was set to zero in the model as the content of dihydrogen in atmosphere is negligible compared to the modeled ones.

The impact of water saturation degree in the soil is not regarded in this study but can have an important impact on hydrogen characteristic transport time in porous media (see e.g. Myagkiy et al., 2020). High water content in the soil reduces its effective porosity and restricts the diffusion and advection of H₂ into the soil. For the shallow subsurface environment up to 30 m down, as in our model, the pressure and temperature are too low to warrant a significant dissolution of hydrogen in water. The dihydrogen arising from deep source is most likely to exceed the limit of its water solubility at these conditions which may lead to formation of the bubbles. That can be one of the ways of transport of dihydrogen through an aquifer. This mode of transport may be described by formation of gas clusters (e.g. Li and Yortsos, 1995; Dominguez et al., 2000) and is out of the scope of this study, which highlights the transport in the vadose zone.

2.2 Boundary conditions

Soil permeability: the average permeability used in most of the simulations reported here is 10 Da (or 10⁻¹¹ m²), appropriate for a moderately permeable soil (Auer et al., 1996; Lazik et al., 2016), the influence of heterogeneity in this anisotropy and especially of lower permeabilities in clayey area will be discussed.

Punctual deep leakage: The modelling case involves a deep origin H₂ source; the gas is, therefore, supposed to vent upward into the overlying rocks and then undergoes gas-soil interactions in the soil layer. The source is represented by a source point at the bottom of the modelling domain, where the gas injects due to overpressure and the dihydrogen enters the porous soil matrix. This periodic concentration pulse of H₂ is obviously an ideal case. However, in the absence of constraints on the emitting zone itself, more complex emission type would be difficult to justify. The Precambrian or ophiolitic rocks, which are the more likely source rock for H₂, being both poorly permeable we will consider that the H₂ emission happens at a given point in the subsurface which corresponds to a fracture that allows, when the fluid pressure is sufficient the escape of the gas. Like a periodic geyser from which gases escape in geothermal area. Other geological contexts may lead to this kind of punctual leakage: the degassing at the top of an aquifer since the H₂ solubility is strongly decreasing with the pressure; leakage by natural hydraulic fracturing is also possible, such a microcrack may get open when local stress tensor changes or at the top of an accumulation if the seal in not good enough as observed above the HC gas accumulation (Sibson, 2003; Petrie, 2014), capillary seal break is another possibility as already proposed for methane leakage (Cartwright and Santamarina, 2015).

Periodicity of the deep H₂ leakage: Since all the available monitoring data show that the H₂ flow is periodic we will mainly use that hypothesis. H₂ monitoring data from Brazil and Kola Peninsula suggest a near 24 h period and a 6 h long flow. We will use these values in our model for the pulsed source. In addition, shorter pulse periods and even continuous H₂ leakage will be modeled for comparison.

Amount of the H₂ leaking: one can suspect that the size of the depression is related to the amount and/or pressure of dihydrogen that is degassing during one cycle. Therefore, the different ΔP_fracture will be tested and the size of the area influenced by the H₂ flow compared to the fairly circles sizes.

Permeability and porosity of the soils: the transport of fluid in soil is highly dependent of its hydraulic characteristics. In the soil, as in the rock, values may be highly variable, of at least three orders of magnitude, from 10 Da (for gravel) to 0.01 Da (to shaly soil). The influence of this parameter and of its distribution within the soil of a same structure will be tested.

Bacterial activity: the impact of this parameter has already been discussed in a previous paper (Myagkiy et al., 2020), so we will incorporate this effect in the discussion but the impact of the variation of the quantity of bacteria and kinetic of the H₂ consumption will not be discussed here in detail.

3 Results and discussion

3.1 Permanent flow of H₂

Figure 2 shows the concentration in the soil above a permanent leakage of H₂ located at 30 m. After a few hours, H₂ reaches the surface and the horizontal concentration profile has a bell shape with a maximum located just above the leak point. Then, the H₂ concentration profile spreads laterally. Therefore, the size of the area affected by the H₂ flow increases with time. Whereas it is about 200 m large after around 2 months, it reaches 400 m after a bit more than a year.

One may note that:

• The H₂ soil concentration has a “bell shape” distribution with its maximum located at the center of the structure

• The requested time for H₂ to reach the surface is around 12 h (number which is obviously depending on the depth of the leakage and pressure of dihydrogen recharge from the deepest level). Above the leak point after 2 days, the [H₂] remains constant in intensity over time. In the center, all the values are so the same, the maximum that could be reached, i.e. the same concentration as at the leakage point

• Progressively, the concentration of the soil is increasing laterally and the affected area, the fairy circle, is getting larger and larger through time. The process of the growing of the structure, is getting slower but continues, the size of the structure, its diameter, is about 50 m after 1 day, 200 m after 50 days and 400 m after 500 days.

Clearly, none of the characteristics of this H₂ signal fits with the observed data described above.

A second set of runs has been performed adding the effect of H₂ consumption in soil due to bacterial activity. The results are shown in Figure 3. In such a case, the kinetics of the consumption bound the growing of the fairy circle. H₂ consumption limits the radial expansion of the H₂ emitting zone. With bacterial activity, after few days the permanent state is achieved, there is an equilibrium between the H₂ emission and consumption, and the size of the structure has reached its maximum which is about 100 m of diameter for this 30-m deep source. The concentration along profile remains thus constant. The balance between the transport time of H₂ in the soil to reach the surface and the kinetics of the H₂ consumption is linked to soil characteristic and bacterial activity; it is absolutely not a one-to-one relationship.

In conclusion, none of these results fit with the literature data, they do not explain the ring of high [H₂] values neither the variabilities of the [H₂] versus time. In addition, one may note that without consumption of H₂ in the soil the distribution of H₂ is continuing and growth of the structure continues over years. In the literature, the soil H₂ concentration is always variable but nobody noticed so far large changes in the structure size through time. Whatever its origin, bacterial or not, H₂ consumption in the soils seems mandatory to explain this stability in case of continuous flow.

Fig. 2 Continuous leakage of H2 with the source at 30 m and without active bacteria in the soil. Soil gas permeability is 10 Da. A. Snapshots of an H2 propagation. The color from blue to red indicates the [H2] of the propagating gas normalized to the source [H2]. B. [H2] distribution as a function of time at 80-cm depth.

Fig. 3 Continuous leakage of H2 with the source at 30 m, considering the active H2-consuming bacteria in the soil. Soil gas permeability is 10 Da. A. Snapshots of an H2 propagation. B. [H2] distribution as a function of time at 80-cm depth.

3.2 H₂ transport reference case pulsed source

A third set of runs has been performed using a pulsed source at a given depth of 30 m, with bacterial activity, the influence of the source depth will be discussed later on. The H₂ pulse chosen as boundary condition at the leakage point is 12 h long and shows a maximum at t = 6 h. Appendix A that resumes the equations and hypotheses shows the concentration of H₂ at the source point for the various runs. Figure 4A shows the rise of the hydrogen from the source, as previously the surface is reached above the leaking point after 12 h. Then, a ring develops defining a high [H₂] area which gets slightly larger until it reaches the maximum diameter of 100 m in the case of a leak point situated at 30 m. Figure 4B shows the same results but highlights the profile versus depth at different times. Interestingly, one may observe that according to the time of the day, the H₂ concentration gradient, $\frac{\partial \left[{\mathrm{H}}_2\right]}{\partial z}$ at the near surface can be negative (t0 + 18 h) or positive (6 h and 12 h) depending on the measurement time.

Figure 4C presents the near surface emanations across the structure at 4 different time-frames, the values are given at 80 cm depth in the center of the structure since it is roughly where the monitoring or the measurement are done in the literature. One may observe the increase of concentration at the center of the structure until 12 h. Then, the peak splits up. The ring forms between 12 and 18 h and then it slightly grows. It is worth to notice that at 24 h (dashed red line), the ring (the two picks in the profile) is perfectly shaped, the [H₂] in the center is zero.

The results of this modeling, even very simple, are a key to understand literature data where apparently random [H₂] values are reported through the structures as summarized in the introduction section. The fact that the source is pulsed induces the ring formation and the inversion of the maximum between the center of the structure and the rim. The delay, i.e. the travel time, between the pulse maximum at the bottom of the studied zone and the moment it reaches the surface explains these versus time and versus depth H₂ concentration curves. The mode of transport, advective and diffusive, and the soil characteristics, essentially its permeability, have a strong influence on the velocity of hydrogen transport as it will be discussed later on.

Fig. 4 Pulse emission of H2 from the source at 30-m depth. A. Snapshots of an H2 gas pulse propagation in time. B. [H2%] distribution along a vertical line passing through the source point (see white line on the t = 0 snapshot of Fig. 4A). Note that the [H2] measured at a given depth will strongly vary over the pulse propagation cycle. C. [H2] distribution at 80-cm depth as a function of time. Reference case: no active bacteria, K = 10 Da, surplus in gas pressure at source is 2 bar, Depth of source = 30 m.

3.3 Size of the structures

The first run with a continuous flow without active bacteria or any other trapping mechanisms already shows that if there is no consumption of H₂ in the soil, H₂ bearing area could be very large, a few hundreds of meters even for a near surface leakage. The scavenging of the H₂ in the soil (e.g., by bacteria) will limit the size of the structure for a given point source flux and depth. The H₂ travel distance which can continuously grow in the absence of H₂ consumption becomes limited when H₂ is consumed on its way to the surface. Accordingly, a deeper source will result in smaller structures or even in no structure at all if all the emitted H₂ is consumed on its way to the surface. However, the literature data in Russia (Larin et al., 2014) as in the USA (Zgonnik et al., 2015) show that in the same area, the size of the structures can be variable. Considering the fact that the thickness of the soil does not usually change too quickly, alternative explanations have to be looked for to explain the variability of the structure diameters. We, therefore, tested the pressure of the H₂ “bubble”, which is prone to vary. Since the time dependency of the pulse is the same (Fig. A1B), that increase of pressure of the H₂ bubble means also than a higher quantity of H₂ is released into the soil from the leak point.

The source being pulsed, the [H₂] signal measured at the surface is time dependent. This time dependency has been described previously so we will focus our analysis on the size of the structure. Therefore, instead of a succession of snapshots, the envelop of oscillating [H₂] will be displayed (Fig. 5), which corresponds to the maximum [H₂] measured on each point of the surface within the first 30 h (with and without bacterial activity).

With the source at 30 m and without bacteria, the circle is 100 m large after 1 day if the pressure is low, 2 bars, and can reach a diameter of 140 m if the H₂ is overpressured (10 bars) when escaping from the basement. These values are reduced to 80 m and 120 m (respectively) in the case of bacterial activity consuming H₂ in the subsurface.

The numerical relationship between size of the structure and depth of the leak point is depending also on the soil permeability and porosity that will be discussed later. However, the prominent observation is the link between the characteristics of the H₂ pulse and the size of the structures. In a given basin, for given soil characteristics and thickness, the size of the structures may vary if the H₂ pulse is more or less pressurized, that means also that the size increases if the quantity of the H₂ released at each pulse is getting larger. It seems particularly interesting to explain the occurence of small structures within the large structures as observed in Carolina bay.

The fact that H₂ consumption in the soil has a strong impact on the size of the structure by limiting its growth has already been shown with the first runs with permanent leakage. The possibility of inverting the surface signal to infer the depth of the leakage point is therefore illusory, too many unknown, or poorly known, characteristics are impacting the size of the emitting zone. In addition to that, we only model the last step of the transport of the gas toward the surface, when the gas invaded the unconsolidated porous media, as evidenced on high resolution seismic data gas chimney (along fault or not) may allow efficient upward migration of gas from deeper level (Gay et al., 2006; Maia et al., 2016). The knowledge of the structural pattern and fault network and so seismic data are mandatory to see where the chimneys root.

Fig. 5 Maximum [H2] envelop at 80-cm depth as a function of source pressure (0.2, 0.5 and 1 MPa) for a 30-meter source (A) without active bacteria and (B) with active bacteria.

3.4 Heterogeneity of the soil

Up to now we have considered a homogeneous soil of constant permeability and porosity in our model; practically it is almost never the case. The circle corresponding to slight depressions, fine particles are transported towards the center by rainwater where they accumulate to form a shaly deposit. It is what has been observed in Russia (Larin et al., 2010) as in Brazil (Prinzhofer et al., 2019). Shaly soils have a very different permeability than silty ones, up to 3 orders of magnitude (e.g. Terzaghi and Peck, 1968; Freeze and Cherry, 1979 in Tab. 2.2, p. 29). To test the influence of permeability heterogeneity on the soil H₂ concentration, we added to our model a low permeability shaly layer at the center, above the emitting point. The computation was performed considering that the low-permeability layer is 3 m thick and 50 m wide and has a permeability of 0.01 Da, i.e. 3 orders of magnitude lower than the surrounding soil, and two times smaller in porosity. In this run the H₂ consumption by bacterial activity is considered.

Figure 6A shows a snapshot of the resulting H₂ flow. In that case, H₂ is first moving upward and laterally as previously but when H₂ flow reaches the low permeability area, it gets around it taking the fastest way to the surface. The ring of high H₂ values already present due to temporal effect is now reinforced at the border of the less permeable area. The maximum is no more at the center of the structure but shifted at the rim of the low permeability layer. This can be clearly seen in Figure 6B that represents the maximum envelop of [H₂] near surface at 80 cm depth, which corresponds to the maximum [H₂] measured at each point of this depth within the first 30 h.

The amplitudes of the signal versus time at the center of the structure and at the border of the shaly layer are shown in Figure 7A, at the same depth as in Figure 6B. It highlights the strong impact of the permeability on the time of H₂ transport. The maximum is reached near the surface on the border after almost 1 day, that means 18 h after the peak of the emission, when it reaches the center after 30 h so 24 h after the peak of the emission. These 6 h of “delay” results in a higher H₂ consumption in the soil and so the signal in the center is close to zero.

It is also worth to note that with bacterial activity 60% of the H₂ is consumed before it reaches the surface but the “travel time” is only slightly affected by the bacteria. The consumption of H₂ affects the H₂ gradient, when it is strongly dependent on the permeability. In that case, after 2 days the signal is stable, the deep H₂ flow reaches the surface after a few hours, up to almost one and a half day for the center or far from the shaly part and then flow and consumption neutralize each other. A steady state is reached.

Fig. 6 A. Snapshots of an H2-bearing gas pulse propagation as for Figure 4 but with the introduction of a low permeability zone (0.01 Da) simulating a clay-rich soil heterogeneity at the center of the structure. Active bacteria are considered here. B. Maximum [H2] envelop at 80-cm depth.

Fig. 7 Concentration at the center and on the border of the structure in case of heterogeneity at the center with K = 0.01 Da and active bacteria. A. Emission every 24 h. B. Emission every 6 h.

3.5 Frequency of the pulse

As already mentioned, a near 24 h long periodicity has been observed within the few published monitoring data, in Russia, in mines, as well as in Brazil in a grass field. The explanation for that period is still debated but the modeling presented here shows that a continuous flow is not compatible with the field data. In addition, we tested the hypothesis of a higher pulse frequency with 4 pulses per days instead of one. In order to compare the results with the previous case all the other parameters are kept the same, the leak point is 30 m deep and there is a low-permeability layer, 50 m large and 3 m thick, in the center of modeled area. Figure 7B shows the [H₂] in the soil in the middle of the structure, so in the shaly part, and at the rim of this shaly part. As in the 24h hours frequency case, Figure 7A, in Figure 7B the H₂ front first reaches the rim, after about 36 h, and then the positions of maximum and minimum are noticed but they are poorly marked. Globally, after 2 days, there is an almost constant signal. It is even more pronounced in the center of the structure, the travel time is larger since the soil permeability is smaller but after 3 days the signal reaches a kind of steady state and the signal is near constant. The amount of H₂ in the soil is globally smaller since the amount of H₂ released at the source is globally smaller. We do not know the exact quantity of H₂ that is leaking in the various basins so we will not discuss so much the final concentrations. However it is worth to realize that even if the soil is just 30 m thick a high frequency pulse will result in almost a permanent signal and oscillation will not be detected by the near surface measurement. After a maximum of 3 days, the H₂ concentration is stable as for a constant source. However, a constant [H₂] at given location in time and with depth does not fit with the observed data.

3.6 Other potential processes that may influence the H₂ concentration in the soil

As explained previously we modeled the transport of H₂ near the surface in an undersaturated high permeability soil. The chosen approach has limitations and doesn’t consider all the physical processes that may exist in the subsurface, but also in the atmosphere, and affect the H₂ flow.

HC gas seeps have been mainly studied offshore and are the marks of the presence of gas in the subsurface; the related pockmarks and/or mud volcano are very common and rather easily imaged with the high-resolution bathymetry and 3D seismic data now available (e.g. Maia et al., 2016 and reference inside). The images clearly show that the gas escape through narrow conduits from the reservoir, which very often, offshore, is the free gas zone below the hydrate layer. These conduits are vertical and may in some cases correspond to small faults even in young, and so uncompacted, sediments. Except onshore near the surface above the piezometric level, all sediments are water-saturated however seismic data as experiments show that the gas may invade unconsolidated porous media as soon as the overpressure overpass the effective stress (Fauria and Rempel, 2011). The hydrogen as the methane is poorly soluble in water at low temperature and low pressure (Lopez-Lazaro et al., 2019). Therefore, we simplify the model and consider only the free gas, not the dissolved one. We also considered that there is enough gas in this area of gas leakage to model it as one single phase and not as isolated bubbles in a water column. Simulation of individual bubbles, their migration and coalescence, could be done at pore scale (e.g. Mahabadi et al., 2018) but at basin scale it is obviously more challenging, and it is not the level of detail that we target in this study. In term of H₂ consumption in the soil, our model is also simple, we only take into account the near surface bacterial activity and the distribution of the bacteria versus depth is an exponential decrease. The H₂ consumption by inorganic processes is not modeled, for instance adsorbtion of H₂ in the shally soil in the center of the structure could also amplify the H₂ signal attenuation phenomenon in the center of the structures but is not computed here.

External causes of H₂ flow variation are also neglected, in our computation the air pressure is set constant over the simulation. Effect of barometric pumping (Auer et al., 1996) is not considered on the advective H₂ transport in the soil. Water content in the soil, due to for instance the seasonality of the rainfall may also change the soil parameters and increase the complexity of the signal.

Despite of all these limitations, of which we are aware, some conclusions can be drawn.

4 Conclusions

The present 2-D model of H₂ transport (advection and diffusion) in a soil layer from a leak point of periodic [H₂] to the surface offers a simple tool to interpret features that are encountered in the field where active emissions of dihydrogen take place. Several parameters may influence the final concentration and the time and depth dependency of this H₂. The goal of this paper is not to adjust a set of parameters to perfectly fit to the data on a given case but rather to understand what controls the 2D distribution and amplitude of hydrogen concentration.

• Observed shapes: fairy circles. Our modeling shows that a punctual hydrogen leakage at a given depth results in a circular structure at the surface. The fact that currently detected emitting structures are not growing and have almost stable shape over few years of observations is incompatible with a permanent flow over years and also strongly suggests the consumption of H₂ in the soil.

• A pulsed source allows to explain the spatial variability of the [H₂] in the soil as observed in Brazil (Prinzhofer et al., 2019) among other localities. The presented model is 2D, but one may note that the fact that the emitting structures are often aligned above accidents as in Russia (Larin et al., 2010) are fully compatible with a model where the H₂ is moving in the basement fractures, escape from the low permeability rock through pumped fracture opening along the main faults and then get transported through another mode discussed here within the soil.

• Temporal variability of the data: the measurements done in the same structure but with a gap of a few days or a few months have shown that the H₂ concentration is not constant and the monitoring done in Brazil have shown that the signal varies all over the day. Moreover, negative gradients such as decrease of signal with depth, and the positive ones where [H₂] does increase are both possible, as was observed in Oman. The modeling has stressed the shift between the [H₂] pulse and the near surface signal even if the soil is perfectly homogeneous. With “regular” hypothesis for the permeability of the soil and the distance to the source point up to few hours of time lag may be noted between the maxima. Few hours between the pulse and the higher concentration near the surface and few hours between the center of the fairy circle and the rim.

• The size of the circle is mainly a function of the transport time between the localized source and the surface. There is so no simple and univocal way to compute the depth of the source point from the size of the circle. For a given set of parameters, the circle will be larger when the source is deeper and has higher recharge pressure but if the permeability of the soil is lower for a given depth the circle will be smaller. Lower permeability leads to bigger contact time with H₂-consuming bacteria. In addition, if a ΔP exists the structure diameter will be also larger.

• Higher values on the rim and near zero in the center of the structure. If the soil is inhomogeneous, the flow of H₂ will be higher in the areas of high permeability. The H₂ flow is slow in the low permeability zone, as a result the gradient is getting smaller and the consumption by bacterial activity is higher. On the field, usually clay minerals fill the center of a depressed structure, a low permeability may form and the H₂ flow may vanish in the center and increase around the low-permeability layer. So pulse and low permeability zone in the center both explain the rim of high values but that the hypothesis of a low permeability on the center that allow to understand why the [H₂] signal remains low in the center which is above the source point.

One may also notice that the so far published data are not compatible with a high frequency pulsed source that lead quickly, after few ten of meters of soil, to a constant signal near the surface.. At this stage, it can be concluded that the inversion of the near surface [H₂] signal in order to deduce the source characteristics (depth, overpressure, [H₂]) is not attainable unless many subsurface parameters are known; among them, soil thickness, permeability variation and H₂ consumption potential are the most important. Permanent monitoring at different depths would also be interesting, for us even mandatory, to go further in the interpretation of the near-surface signal.

Acknowledgement

We thank ENGIE who supported this study and Andrey Myagkiy post-doctoral position. Complementary funding by CNRS/INSU through the TelluS Program is also acknowledged.

Appendix A Physical governing equations

Equations that govern the hydrodynamic system are represented by the equation of gas transport (1), continuity equation (2) and Darcy flow equation (3). Darcy’s Law is combined with continuity equation and the ideal gas equation (4) becomes equation (5). $$\frac{\partial \left(\phi c_i\right)}{\partial t}+\nabla \cdot \left(u{c}_i\right)=\nabla \cdot \left(D^{*}\nabla c_i\right)+R^{*}\text{,}$$(1) $$\frac{\partial \left(\rho \phi \right)}{\partial t}=-\nabla \cdot \left(\rho -\frac{k}{\mu }\nabla P\right)\text{,}$$(2) $$u=-\frac{k}{\mu }\left(\nabla P-\rho g\right)\text{,}$$(3) $$\rho =\frac{PM}{RT}\text{,}$$(4) $$\frac{\partial \left(P\phi \right)}{\partial t}=-\nabla \cdot \left(-\frac{k}{\mu }P\nabla P\right)\text{,}$$(5) $$D_{por}=\phi \tau D_{atm,}$$(6) $$\tau ={\phi }^{1/3}\text{,}$$(7) $$D^{*}=D_{por}+\alpha \left|u\right|\text{.}$$(8)

For all-gas conditions as discussed in this paper, the tortuosity factor reduces to the porous medium value (Eq. (7)), according to the tortuosity correlation of Millington and Quirk (1961).

Therefore, the system of equations is the following:

where, c_i is the concentration of H₂ in a gas phase within porous media, R^* is a rate of H₂ consumption by bacteria, and u signifies the gas velocity, ϕ denotes the porosity, τ is the tortuosity, D^* is the sum of effective diffusion coefficient for the hydrogen gas in a porous media and a velocity-dependent dispersion term, so called mechanical dispersion, and D_a is a diffusion coefficient of H₂ in air.

The diffusive transport in the gas was traditionally described by molecular diffusion. The diffusive flux was modeled with the use of Fick’s law that has been optimized for porous media and takes into account porosity and tortuosity factors. The molecular mean free path of gas particles was assumed to be way larger than pore dimension and, thus, no influence of pore walls was considered. Preferential transport pathways in the soil, such as presence of cracks, worm channels or decayed and live root channels have not been considered in this model.

The modelling case involves a deep origin H₂ source; the gas is, therefore, supposed to vent upward into the overlying sediments and then undergo the gas-soil interactions in the soil horizon. The source is represented by a source point at the bottom of the modelling domain, where the gas injects due to the pressure surplus and the dihydrogen enters the porous matrix of domain, the boundary conditions versus time is presented in Figure A1.

Fig. A1 Signal at the source leaking point. A. Permanent source. B. Periodic source with one emission a day. C. Periodic source with four emissions a day.

For numerical calculations reported here, the permeability used is 10 Da (or 10⁻¹¹ m²), appropriate for a moderately permeable soil. Lower permeabilities should cause smaller effects.

The main physical parameters used in the model are listed in Table A1.

References

• Abrajano TA, Sturchio NC, Bohlke JK, Lyon GL, Poreda RJ, Stevens CM. 1988. Methane-hydrogen gas seeps, Zambales Ophiolite, Philippines: deep or shallow origin? Chemical Geology 71(1-3): 211–22. DOI: 10.1016/0009-2541(88)90116-7.
• Auer LH, Rosenberg ND, Birdsell KH, Whitney EM. 1996. The effects of barometric pumping on contaminant transport. Journal of Contaminant Hydrology 24(2): 145–166.
• Brunet F. 2019. Hydrothermal production of H₂ and magnetite from steel slags: a geo-inspired approach based on olivine serpentinization. Front. Earth Sci. 7. DOI: 10.3389/feart.2019.00017.
• Cartwright J, Santamarina C. 2015. Seismic characteristics of fluid escape pipes in sedimentary basins: implications for pipe genesis. Mar. Petroleum Geol. 65: 126e140.
• Charlou JL, Donval JP, Fouquet Y, Jean-Baptiste P, Holm N. 2002. Geochemistry of high H₂ and CH₄ vent fluids issuing from ultramafic rocks at the Rainbow hydrothermal field (36°14′N, MAR). Chemical Geology 191(4): 345–59.
• Charlou J-L, Donval JP, Konn C, Ondreas H, Fouquet Y, Jean-Baptiste P, et al. 2010. High production and fluxes of H₂ and CH₄ and evidence of abiotic hydrocarbon synthesis by serpentinization in ultramafic-hosted hydrothermal systems on the mid-atlantic ridge. Geophysical Monograph Series 188. Washington DC: American Geophysical Union.
• Conrad R, Seiler W. 1980. Contribution of hydrogen production by biological nitrogen fixation to the global hydrogen budget. Journal of Geophysical Research Oceans 85: 5493–5498.
• Deronzier JF, Giouse H. 2020. Vaux-en-Bugey (Ain, France): the first gas field produced in France, providing learning lessons for natural hydrogen in the sub-surface? BSGF – Earth Sci. Bull. 5: 191. DOI: 10.1051/bsgf/2020005.
• Deville E, Prinzhofer A. 2016. The origin of N₂–H₂–CH₄-rich natural gas seepages in ophiolitic context: A major and noble gases study of fluid seepages in New Caledonia. Chem. Geol. 440: 139–47.
• Dominguez A, Bories S, Prat M. 2000. Gas cluster growth by solute diffusion in porous media. Experiments and automaton simulation on pore network. International Journal of Multiphase Flow 26(12): 1951–1979.
• Fauria K, Rempel A. 2011. Gas invasion into water-saturated, unconsolidated porous media: implication for gas hydrate reservoirs. EPSL 312: 188–193. DOI: 10.1016/j.epsl.2011.09.042.
• Firstov PP, Shirokov VA. 2005. Dynamics of molecular hydrogen and its relation to deformational processes at the Petropavlovsk-Kamchatskii geodynamic test site: evidence from observations in 1999–2003. Geochemistry International. 43(11): 1056–1064. [Translated from Geokhimiya 2005(11): 1151–1160].
• Freeze RA, Cherry JA. 1979. Groundwater. NJ: Prentice Hall, 604 p.
• Gay A, Lopez M, Cochonat P, Levaché D, Sermondadaz G, Seranne M. 2006. Evidences of early to late fluid migration from an Upper Miocene turbiditic channel revealed by 3D seismic coupled to geochemical sampling within seafloor pockmarks, Lower Congo Basin. V 23(3): 387–399. DOI: 10.1016/j.marpetgeo.2006.02.004.
• Guélard J, Beaumont V, Rouchon V, Guyot F, Pillot D, et al. 2017. Natural H₂ in Kansas. Deep or shallow origin? Geochem. Geophys. Geosyst. 18(5): 1841–65.
• Klein F, Tarnas J, Bach W. 2020. Abiotic sources of molecular hydrogen on Earth. Elements 16: 19–24.
• Larin NV, Larin VN, Gorbatikov AV. 2010. Circular structures, caused by the deep seeping of hydrogen. Degassing of the Earth: Geotectonics, geodynamics, Deep fluids, oil and gas, hydrocarbons and life. Russian conference with international participation dedicated to the 100th anniversary of academician P.N. Kropotkin, 282 p. (in Russian)
• Larin N, Zgonnik V, Rodina S, Deville E, Prinzhofer A, Larin VN. 2014. Natural molecular hydrogen seepage associated with surficial, rounded depressions on the European Craton in Russia. Natural Resources Research 24(3): 369–83.
• Lazik D, Ebert S, Neumann PP, Bartholmai M. 2016. Characteristic length measurement of a subsurface gas anomaly – Amonitoring approach for heterogeneous flow path distributions. International Journal of Greenhouse Gas Control 47: 330–341.
• Li X, Yortsos YC. 1995. Theory of multiple bubble growth in porous media by solute diffusion. Chemical Engineering Science 50(8): 1247–1271.
• Lopez-Lazaro C, Bachau P, Moretti I, Ferrando N. 2019. Hydrogen solubility in aqueous NaCl solutions: from molecular simulation to equation of state. BSGF. DOI: 10.1051/bsgf/2019008.
• Mahabadi N, Zheng X, Yun TS, van Paassen L, Jang J. 2018. Gas bubble migration and trapping in porous media: Pore-scale simulation. Journal of Geophysical Research: Solid Earth 123: 1060–1071. DOI: 10.1002/2017JB015331.
• Maia A, Cartwright J, Andersen E. 2016. Shallow plumbing system inferred from spatial analysis of pockmark array. Marine and Petroleum Geology 77: 865–881. DOI: 10.1016/j.marpetgeo.2016.07.029.
• Mayhew LE, Ellison ET, McCollom TM, Trainor TP, Templeton AS. 2013. Hydrogen generation from low-temperature water–rock reactions. Nature Geoscience 6: 478–484. DOI: 10.1038/ngeo1825.
• Miller HM, Mayhew LE, Ellison ET, Kelemen P, Kubo M, Templeton AS. 2017. Low temperature hydrogen production during experimental hydration of partially-serpentinized dunite. Geochimica et Cosmochimica Acta 209: 161–183. DOI: 10.1016/j.gca.2017.04.022.
• Millington RJ, Quirk JP. 1961. Permeability of porous solids. Transactions of the Faraday Society 57: 1200–1207.
• Monnin C, Chavagnac V, Boulart C, Ménez B, Gérard M, Gérard E, et al. 2014. Fluid chemistry of the low temperature hyperalkaline hydrothermal system of Prony Bay (New Caledonia). Biogeosciences 11(20): 5687–5706. DOI: 10.5194/bg-11-5687-2014.
• Myagkiy A, Brunet F, Popov C, Krüger R, Guimarães H, Charlet L, et al. 2020. H₂ dynamics in the soil of an H2-emitting zone (São Francisco Basin, Brazil): Microbial uptake quantification and reactive transport modelling. Applied Geochemistry 212. DOI: 10.1016/j.apgeochem.2019.104474.
• Neal C, Stanger G. 1983. Hydrogen generation from mantle source rocks in Oman. Earth and Planetary Science Letters 66: 315–320.
• Nivin VA, Pukha VV, Lovchikov AV, Rakhimov RG. 2016. Changes in the molecular hydrogen concentration in an underground mine (Lovozero rare-metal deposit, Kola peninsula). Doklady Earth Sciences 471(2): 1261–1264.
• Noel J, Godard M, Oliot E, I Martinez. 2018. Evidence of polygenetic carbon trapping in the Oman Ophiolite: Petro-structural, geochemical, and carbon and oxygen isotope study of the Wadi Dima harzburgite-hosted carbonates (Wadi Tayin massif, Sultanate of Oman). Lithos 323: 218–237.
• Petrie E. 2014. Rock strength of caprock seal lithologies: Evidence for past seal failure, migration of fluids and the analysis of the reservoir seal interface in outcrop and the subsurface. Dissertation Thesis, Utah State University, Logan, Utah, 208 p.
• Prinzhofer A, Deville E. 2015. Hydrogene naturel. La prochaine révolution énergétique ? Paris : Belin.
• Prinzhofer A, Cisse CST, Diallo AB. 2018. Discovery of a large accumulation of natural hydrogen in Bourakebougou (Mali). Int. J. Hydrogen Energy 43: 19315–26.
• Prinzhofer A, Moretti I, Françolin J, Pacheco C, D’Agostino A, Werly J, et al. 2019. Natural hydrogen continuous emission from sedimentary basins: The example of a Brazilian H2-emitting structure. International Journal of Hydrogen Energy 44(12): 5676–5685. DOI: 10.1016/j.ijhydene.2019.01.119.
• Sano Y, Urabe A, Wakita H, Wushiki H. 1993. Origin of hydrogen–nitrogen gas seeps, Oman. Applied Geochemistry 8: 1–8.
• Sato M, Sutton AJ, McGee KA, Russel-Robinson S. 1986. Monitoring of hydrogen along the San Andreas and Calaveras Faults in Central California in 1980–1984. Journal of Geophysical Research 91(B12): 12315–12326.
• Semlitsch RD. 2000. Size does matter: the value of small isolated wetlands. Natl Wetl Newsl 22: 5–6.
• Sherwood Lollar B, Onstott TC, Lacrampe-Couloume G, Ballentine CJ. 2014. The contribution of the Precambrian continental lithosphere to globa H₂ production. Nature 516: 379–382. DOI: 10.1038/nature14017.
• Sibson RH. 2003. Brittle-fault controls on maximum sustainable overpressure in different tectonic regimes. AAPG Bulletin 87: 901–908. DOI: 10.1306/01290300181.
• Terzaghi K, Peck RB. 1968. Soil mechanics in engineering practice. New York: John Wiley & Sons, 729 p.
• Vacquand C, Deville E, Beaumont V, Guyot F, Sissmann O, Pillot D, et al. 2018. Reduced gas seepages in ophiolitic complexes: Evidences for multiple origins of the H₂–CH₄–N₂ gas mixtures. Geochimica Cosmochimica Acta 223: 437–461.
• Welhan JA, Craig H. 1979. Methane and hydrogen in East Pacific rise hydrothermal fluids. Geophys. Res. Lett. 6: 829–831.
• Zgonnik V. 2020. The occurrence and geoscience of natural hydrogen: A comprehensive review. Earth Science Review 203. DOI: 10.1016/j.earscirev.2020.103140.
• Zgonnik V, Beaumont V, Deville E, Larin N, Pillot D, Farrell K. 2015. Evidences for natural hydrogen seepages associated with rounded subsident structures: the Carolina bays (Northern Carolina, USA). Prog. Earth Planet. Sci. 2: 31. DOI: 10.1186/s40645-015-0062-5.

All Tables

All Figures

Fig. 1 Synthesis of the fairy circle emissions characteristics. A. A fairy circle in Brazil showing how the vegetation is affected. B. Variability on time over 5 weeks for a given sensor in the structure monitored in Brazil. See Prinzhofer et al. (2019) and Myagkiy et al. (2020) for details about that case. C. Daily variations for 6 sensors installed on the structure monitored in Brazil showing the 24 h period. D. Variability in size, of the emitting structures in USA (modified from Prinzhofer and Deville, 2015). E. Two profiles across a structure showing the spatial variability of H2 concentration value between the center and the rim of the structure, Lake Podowoye, Russia. These values have not been taken at the same time but over 2 days. Modified from Larin et al. (2014). F. Profile showing an example of [H2] concentration variation versus depth in USA (modified from Zgonnik et al., 2015).

In the text

	Fig. 2 Continuous leakage of H2 with the source at 30 m and without active bacteria in the soil. Soil gas permeability is 10 Da. A. Snapshots of an H2 propagation. The color from blue to red indicates the [H2] of the propagating gas normalized to the source [H2]. B. [H2] distribution as a function of time at 80-cm depth.
In the text

	Fig. 3 Continuous leakage of H2 with the source at 30 m, considering the active H2-consuming bacteria in the soil. Soil gas permeability is 10 Da. A. Snapshots of an H2 propagation. B. [H2] distribution as a function of time at 80-cm depth.
In the text

Fig. 4 Pulse emission of H2 from the source at 30-m depth. A. Snapshots of an H2 gas pulse propagation in time. B. [H2%] distribution along a vertical line passing through the source point (see white line on the t = 0 snapshot of Fig. 4A). Note that the [H2] measured at a given depth will strongly vary over the pulse propagation cycle. C. [H2] distribution at 80-cm depth as a function of time. Reference case: no active bacteria, K = 10 Da, surplus in gas pressure at source is 2 bar, Depth of source = 30 m.

In the text

	Fig. 5 Maximum [H2] envelop at 80-cm depth as a function of source pressure (0.2, 0.5 and 1 MPa) for a 30-meter source (A) without active bacteria and (B) with active bacteria.
In the text

	Fig. 6 A. Snapshots of an H2-bearing gas pulse propagation as for Figure 4 but with the introduction of a low permeability zone (0.01 Da) simulating a clay-rich soil heterogeneity at the center of the structure. Active bacteria are considered here. B. Maximum [H2] envelop at 80-cm depth.
In the text

	Fig. 7 Concentration at the center and on the border of the structure in case of heterogeneity at the center with K = 0.01 Da and active bacteria. A. Emission every 24 h. B. Emission every 6 h.
In the text

	Fig. A1 Signal at the source leaking point. A. Permanent source. B. Periodic source with one emission a day. C. Periodic source with four emissions a day.
In the text