Deposition velocity of H2: a new algorithm for its dependence on soil moisture and temperature

1. Introduction

The major sink of atmospheric molecular hydrogen, H₂, is the uptake by soil. This uptake totals about 60 Tg H₂/yr and is a widespread process: Virtually, all soils containing organic carbon were shown to take up H₂. Moreover, the uptake was shown to be predominantly biological – mediated by enzymes, i.e. by various kinds of hydrogenase (see Conrad, 1999; Constant et al., 2010), and to take place in the uppermost 10 cm of the soil (Liebl and Seiler, 1976; Förstel, 1986; Yonemura et al., 1999).

The rate of uptake is usually expressed as deposition velocity, v_d [cm/s], such that the flux of H₂ from the atmosphere to the soil surface, F_a [molec/cm³/s], is given by

1

This deposition velocity was measured by numerous in situ experiments for various biomes and shown to vary between 0.01 and 0.1 [cm/s] (see Ehhalt and Rohrer, 2009, for a summary). The field measurements also demonstrated a significant decrease of v_d with increasing soil moisture except occasional observations of very low v_d for very dry soils. The temperature dependence of v_d measured in the field was relatively weak and exhibited a broad maximum around 30°C (Liebl and Seiler, 1976; Conrad and Seiler, 1985; Förstel, 1986; Yonemura et al., 1999; Yonemura et al., 2000a, 2000b; Lallo et al., 2008; Schmitt et al., 2009). Unfortunately, the measurements so far do not suffice to characterize all relevant types of soil, and the estimates of global average uptake based on them involve large extrapolations.

More recently, there have also been several attempts to formulate theoretical expressions for v_d (or of F_a directly), which allow us to predict its dependence on soil moisture and temperature (Yonemura et al., 2000b; Smith-Downey et al., 2008; Schmitt et al., 2009; Bousquet et al., 2011). They were based on 1D molecular diffusion of H₂ into and within the soil and on its enzymatic removal within. To achieve simple expressions, it was usually assumed that the soil parameters (soil porosity, moisture, temperature, as well as enzymatic activity) are uniform throughout the soil. Thus, most of these studies disregarded the observation that the top layer of the soil is often so dry that its enzymatic activity is greatly diminished or even totally destroyed (Yonemura et al., 2000b; Smith-Downey et al., 2008). Recognizing this, Yonemura et al. (2000b) therefore proposed a two-layer soil diffusion model for the uptake of H₂. The upper layer is assumed to be devoid of hydrogenase activity and acts essentially as a diffusion barrier. Its thickness, δ, on the order of 1 cm, increases with dryness. In the layer below, both diffusion and removal of H₂ take place. The authors solved the model numerically, and the results matched reasonably well the observed dependence of v_d on soil moisture.

Here we will show that using this model and virtually the same assumptions about soil uniformity as Yonemura et al. (2000b), we can obtain an analytical solution for v_d (Section 2). This solution explicitly depends on δ.

To characterize the enzymatic removal rate of H₂ within the soil and its dependence on soil moisture and temperature, we rely on published data from laboratory experiments (Section 3). With an expression for the average dependence of δ on soil moisture, which is derived in Section 4, we can formulate a complete dependence of v_d on soil moisture and temperature. In Section 5, the v_d calculated from our analytical expression is compared against those disregarding the effects of a dry layer and a set of v_d from field measurements. The implications are discussed.

2. Derivation of a closed expression for v_d

2.1. Diffusion equations

Following the earlier studies (Yonemura et al., 2000b; Smith-Downey et al., 2008), we assume that molecular diffusion drives the transport of H₂ within the soil. We further assume that there is no production of H₂ in the soil and that the destruction is first order in the H₂ concentration. Therefore, the 1D diffusion equations for flux and mass balance take the form:

2

3

where F_s is the flux in molec cm⁻² s⁻¹, ρ is the number density of soil air in molec cm⁻³, M_s is the mixing ratio of H₂ in soil air, and z is the depth in cm. k_s is the rate constant for removal of H₂ from soil air in s⁻¹. k_s ⁻¹=τ, the mean lifetime of H₂ in soil air.

Θ_a is the fraction of soil volume filled with air. It is given by

4

where Θ_p is the total porosity of the soil (cm³ total pore space/cm³ soil), and Θ_w is the volumetric soil water content (cm³ of water-filled pores/cm³ soil).

D_S is the diffusivity of H₂ in the soil (units: cm² s⁻¹). There are several ways to formulate D_S. Following Millington and Quirk (1961) (see also Yonemura et al., 2000b; Schmitt et al., 2009), D_S can be expressed as:

5

where D_A is the molecular diffusivity of H₂ in air (units: cm² s⁻¹). D_A depends on temperature, T in °C, and atmospheric pressure, p in hPa (see Yonemura et al., 2000b):

6

Our calculations will assume steady state, i.e.

7

2.2. Two-layer model

To calculate the deposition velocity, we apply these diffusion equations to a two-layer soil model first proposed by Yonemura et al., 2000b. It consists of an upper dry layer (I) of depth δ without H₂ removal and a more moist layer (II) below, in which H₂ is destroyed enzymatically at a rate characterized by the uniform rate constant k_s. Since layers I and II differ in soil moisture, Θ_w, they also differ in diffusivity, D_S [see eq. (5)]. The deposition velocity is defined by eq. (1) given above to

8

In layer l the hydrogen flux, F_I, remains constant and equal to F_a. In that layer, i.e. for 0≤z≤δ, the gradient in M_s is therefore linear and given by

9

where M_s(0)=M_a and M_s(δ) are the H₂ mixing ratios in the soil at depths z=0, and z=δ, respectively.

For z≥δ, i.e. layer II, the vertical profile of M_s is determined by eq. (7). Assuming D_SII, ρ, Θ_a, k_s to be constant with depth, eq. (7) can be solved analytically:

10

The H₂ flux at depth δ can be calculated by inserting solution (10) in eq. (2)

11

Or solving for M_s(δ)

12

Inserting eqs. (12) and (8) into eq. (9), and observing that F_a=F_I=F_II(δ), and M_a=M_s(0), we solve for v_d and obtain

13

Eq. (13) presents an analytical expression for v_d. Besides the simple soil parameters, D_SI, D_SII, Θ_a, it also contains the more involved quantities δ and k_s. The latter also depend on the soil moisture; k_s also on the soil temperature. Thus, before eq. (13) can be used to describe the dependence of v_d on soil moisture and temperature, expressions for those individual dependencies must be found. This will be done in the two following sections.

In passing, we note that eq. (13) unifies two special cases that have been used earlier to parameterise v_d. For moist soils, δ approaches 0, and eq. (13) reduces to the solution as derived by Yonemura et al. (2000b) for their one-layer model (see also Smith-Downey et al., 2008). For very large δ, i.e. , eq. (13) reduces to V_d=D_SI/δ, the expression used by Schmitt et al. (2009).

3. The dependence of k_sΘ_a on soil temperature and soil moisture

k_s and its dependence on soil moisture and temperature can and has been determined by laboratory experiments. To be exact, what these experiments have determined is the functional form of the dependence of k_sΘ_a on temperature and moisture (see Ehhalt and Rohrer, 2011). Moreover, as the experiments of Smith-Downey et al. (2006) demonstrated, the dependences on soil moisture and temperature can be separated, i.e., k_sΘ_a can be written as

14

15

For the moisture dependence f_es(Θ_w/Θ_p) of eolian sand and f_ll(Θ_w/ Θ_p) for loess loam they found

16

17

The only other soil for which f(Θ_w/Θ_p) has been published is boreal forest soil (Smith-Downey et al., 2006), which has also been reanalysed by Ehhalt and Rohrer (2011). The function g(T) is depicted in Fig. 1; the functions f_es(Θ_w/Θ_p) and f_ll(Θ_w/Θ_p) are shown in Fig. 2. For more details, see Ehhalt and Rohrer (2011).

Fig. 1. 

Fig. 2. 

4. The dependence of δ on soil moisture

The depth δ is an auxiliary quantity. It captures with a single parameter the essence of the influence of the vertical soil moisture profile on v_d. This simple parameterisation seems justified, since well-resolved vertical profiles of soil moisture are hardly ever reported when v_d is measured. Under the best circumstances, the average soil moisture across the top 10 cm of soil is available – usually measured by time domain reflectometry. Thus, to enable comparison with field experiments, we choose this average, , as a measure of soil moisture.

To derive an expression for the dependence of δ on , we proceed semi-empirically. We first calculate sequences of realistic soil moisture profiles with the help of the 1D model of HYDRUS.¹ [1] The calculations assume a 4 m deep uniform soil column with a total porosity, Θ_p=0.38, and temperature, T=15°C, bounded at the lower end by the water table. The vertical resolution dh varies from 0.008 cm at the upper bound to 0.8 cm at the lower bound according to (dh)²=h/625 cm⁻¹+0.008² cm². For realistic upper boundary conditions, we select the measured daily mean values of rainfall, air temperature, and relative humidity obtained at the weather station of the Jülich Research Centre during the year 2008. Two types of soil are investigated: Sandy loam, and loam. For technical reasons, the profiles are read out every 1.5 days.

From eq. (14) and the calculated soil moisture profiles, we generate vertical profiles of the removal rate, k_sΘ_a, from which we will eventually estimate δ.

To illustrate the types of profiles obtained, Fig. 3 describes the sequences of profiles of soil moisture and removal rate in sandy loam and loam during a dry spell of about two weeks following a day with rainfall (27 April 2008). As the sequence for sandy loam indicates, the soil moisture profile develops a shallow dry layer with Θ_W close to zero and about 0.1 cm depth within 1.5 days. With increasing time, this layer deepens and reaches about 1 cm depth after two weeks. The dry layer maintains a finite level of soil moisture which depends somewhat on the model assumptions. The exact value of that moisture is, however, not important, as long as it remains lower than the threshold moisture below which the enzymatic H₂ uptake vanishes. This dry layer is bounded below by a narrow zone of rapidly increasing moisture followed by a more gradual increase in moisture at greater depth. The zone of rapid increase in soil moisture has been termed ‘evaporation front’ or ‘drying front’ to stress the fact that soil water transport proceeds in the form of water vapour above and in the form of liquid water below the front (e.g. Saravanapavan and Salvucci, 2000).

Fig. 3. 

The corresponding profiles of the removal rate constant, k_sΘ_a, are similar to those for soil moisture. For sandy loam, the rapid increase is even accentuated owing to the fact that the dependence of k_sΘ_a on Θ_w for sand, eq. (16) shows a steep rise at low Θ_w (see Fig. 2). As a consequence δ is quite well characterized. In the following, δ will be defined as the depth at which the ½ – value of the maximum k_sΘ_a in the top 10 cm of soil is reached. Panel (3a) also indicates that k_sΘ_a in a given profile is fairly uniform at greater depths. Hence, it corroborates the two-layer model which approximates the actual k_sΘ_a profiles by a Heaviside function with a jump of k_sΘ_a from 0 to a finite, constant value at depth δ.

A similar picture emerges for loam – except that the soil remains moister and the dry zone shallower (Fig. 3b). In addition, the depth profiles for k_sΘ_a are more gradual owing to the slower increase of k_sΘ_a with soil moisture for loam, eq. (17).

The δ so obtained from the various annual series of profiles are plotted against the corresponding in Fig. 4, where is obtained from integrating the soil moisture profiles over the first 10 cm. The total ensemble of data points follows similar hyperbolical curves. However, Fig. 4 also indicates that the individual soils actually have their own hyperbolical shapes with the curves for the soils of smaller pore size ranging over higher average soil moistures, . The full curves represent least square fits to the data for sandy loam, eq. (18), and loam, eq. (19):

18

19

Fig. 4. 

Eq. (18) should hold for all sandy soils. In fact, we also calculated δ for pure sand, but the data clustered too closely for a meaningful fit. Nevertheless, they were compatible with eq. (18). Moreover, calculations for loamy sand were reasonably close to that for sandy loam.

In the choice of the functions and their argument, (Θ_p – )/, we were guided by the observations that: (1) δ=0 when all the soil pores are filled with water, i.e. when =Θ_p; (2) δ approaches 10 cm when approaches the level , the soil moisture where the H₂ removal rate vanishes; (3) δ() is monotonic. The latter suggests a power function or a power series for δ(). For simplicity we choose a power function. Obviously, the derived δ() depends on the choice of Θ_p and k_s(T, Θ_w) Θ_a, but mostly it depends on the type of soil chosen in modelling the depth profile of Θ_W.

We also note that the actual data points scatter substantially around the average curves for δ_S() and δ_L(). The reason is, of course, that different profiles of soil moisture – leading to different δ - can result in the same . In particular, δ-values of zero are reached even at the lowest . This occurs when soil moisture profiles are sampled during the time of rainfall. Then the soil moisture exceeds even in the topmost soil layer. Therefore, δ depends not only on but also on the immediate history of rainfall and evaporation. It also means that the use of eqs. (18) and (19) will introduce an additional uncertainty in the eventual calculation of .

5. Discussion

With the help of eq. (14) for k_sΘ_a, and eq. (18) or eq. (19) for δ, we can calculate v_d from eq. (13) for sandy or loamy soils, if Θ_p, and T are known. In particular, we can calculate the dependence of v_d on for such soils. Figure 5 presents examples for sandy soil and loam, and compares them with the corresponding v_ds from the single-layer model, where . For the later comparison with field measurements of v_d, the shown dependences are based on Θ_p=0.38, T=15°C. Moreover, the calculations assume eq. (16) for the moisture dependence of k_sΘ_a in sand and eq. (17) for that in loam. For both soils, the parameter A in eq. (14) was set to 1, and the soil water contents Θ_WI in the top layer δ to , where is the threshold moisture below which k_sΘ_a vanishes. For Θ_p=0.38, it is 0.01 for sandy soil, and 0.02 for loam. Θ_WII is then calculated by placing the remaining, major fraction of in the soil layer between δ and 10 cm:

20

Fig. 5. 

Figure 5 demonstrates three major features in the dependence of v_d on . (1) All profiles exhibit the same general pattern, namely v_d=0 for , v_d approaching 0 when approaches Θ_p, and a more or less broad maximum in between. (2) The form of v_d() differs significantly for different soils. This is caused mainly by the difference in k_sΘ_a(Θ_w). (3) The inclusion of a dry top layer introduces a significant decrease in the calculated v_d at smaller . That decrease can amount to up to 40% and depends on the type of soil.

The validity of v_d predicted from eq. (13) can be tested: (1) By comparison with existing field measurements of v_d in soils, for which , Θ_p, and T have also been measured. (2) By comparison with the numerical results from modelling an idealized soil. In the second case not only Θ_p and T can be specified, but also vertical profiles of Θ_w and k_sΘ_a in the soil, providing a much more detailed characterization of the soil conditions. In fact, with such information v_d can be calculated directly by numerically solving the differential eq. (7). With the data on vertical profiles of Θ_w and k_sΘ_a acquired for deriving δ in Section 4 we can exactly do that. Figure 6 shows the results for sandy loam and loam. The full curves are the two-layer model cases from Fig. 5 and represent the average dependence of v_d on . They agree well with the data points calculated numerically from the individual profiles of Θ_w and k_sΘ_a. The agreement with the one-layer model is significantly less. This test clearly favours the two-layer model.

Fig. 6. 

However, we note that even in that case, the individual data points can deviate substantially from the mean, a fact that was already observed for δ (Fig. 4) and assigned to the observation that quite different profiles of Θ_w can lead to the same . We further note that for v_d the deviations tend to become larger for larger . A likely explanation is that longer periods without rainfall are required to generate dry soils. During that time the moisture profiles reach similar and regular shapes (see Fig. 3). In contrast, during periods of frequent rainfall events and therefore higher water content, the moisture profiles maintain more irregular and varied shapes.

The deviations of the individual v_d from the mean predicted by eq. (13) (Fig. 6) can be used to estimate the possible fluctuations of individual v_d measured in the field caused by variability that is unaccounted for in the moisture profiles. The relative deviations for loam, δ_vd, show an average moisture dependence well approximated by the relation:

21

Fig. 7. 

Clearly, the v_d measured in the field reach values more than a factor of two higher than those calculated for loam in Fig. 6b. This indicates that the uptake rate constant, k_s(Θ_w,T)Θ_a, in the bare loess loam was higher than that assumed in Fig. 6b, and suggests that the scaling factor A in eq. (14) must be significantly larger than 1 to fit the measured data. We further note that the scatter in the experimental v_d appears even larger than that observed in Fig. 6b. One explanation is that k_s(Θ_w,T)Θ_a varied with time. This is most likely because the temperature varied with season. A second possibility is that the enzymatic activity varied from site to site. This possibility is also quite plausible considering that the v_d from site 1 — indicated by the open triangles in Fig. 7 – all tend to fall into the upper range of the measured v_d().

To accommodate the latter possibility, we present two choices of v_d() calculated from eq. (13): One with A=2, the other with A=8. The heavy line assumes A=2, eqs. (14), (15), and (17) to characterize k_sΘ_a(,T), eq. (19) for δ, Θ_p=0.38, and T=15°C. To provide a sense of the spread in v_d introduced by the range in observed temperatures or by the variance of the soil moisture profiles leading to the same , it is accompanied by error bounds. The dark area gives the range of v_d due to the range in T. It is obtained by calculating v_d from eq. (13) for T=0°C and 30°C, respectively, but otherwise the conditions are identical. (The relatively small spread with temperature apparent in the calculated v_d also offers an explanation why the v_d observed in the field show a large dependence on Θ_w but a weak one on T.) The light shaded area covers the deviations due to the variance in the soil moisture profiles. It is obtained by multiplying v_d() for 30°C with () from eq. (21) and by dividing v_d() for 0°C by that term.

The thin line makes the same assumptions as that for A=2, except that A=8 and T=30°C; the latter to account for the fact that the largest v_d was obtained at temperatures≥30°C.

Obviously, the comparison of the modelled v_d() with this set of field data is not conclusive. On the one hand, the upper and lower most of the measured v_d appear to follow envelopes which are compatible with the modelled dependence of v_d on . In particular, v_d() for A=2, and its error bounds cover the field data – apart from those from site 1 – reasonably well. This is gratifying, since apart from the factor A, the calculated v_d() is not based on the field data for v_d, but are derived from independent information: Laboratory measurements for k_s(Θ_w,T)Θ_a and model calculations for δ.

On the other hand, the field data on v_d as a whole do not show the moisture dependence predicted by eq. (13). In fact, their dependence is most accurately fitted by a straight line. For the data without those from site 1, the regression line takes the form v_d=0.065 – 0.2*, with a standard deviation of 0,009 cm/s. For comparison, the standard deviation of these data from the v_d() predicted from eq. (13) with A=2 is 0.013—clearly a worse fit.

There are several reasons for this: (1) The present field data do not reach to low enough to produce an unequivocal turnaround in the measured v_d. (2) The data are not homogeneous enough. They are possibly influenced by a variation in k_s(Θ_w,T)Θ_a from site to site, which makes it impossible to describe them with one single factor A or a single temperature. (3) They are not numerous enough to select sufficiently large subsets of v_d measured at the same site and at the same temperature and thus to remove some of the variance in the v_d. (4) Eq. (17) for k_sΘ_a() used to model v_d might not exactly match the k_sΘ_a() acting in the measured soils. Indeed a superposition of eqs. (16) and (17) corresponding to a small admixture of sand would shift the maximum in the resulting k_sΘ_a() to lower and yield a better fit of v_d() from eq. (13) to the field data.

For reasons 1–3, these field data also do not lend themselves to differentiate between the v_d() predicted from a one-layer and a two-layer model.

The differences between v_d() from a one-layer and a two-layer model vary not only with the type of soil (Fig. 6), but also with the enzymatic activity in the same type of soil. The latter is demonstrated in Fig. 8 which compares the ratios of v_d(two layer) to v_d(one layer) as a function of for loam with various values of A. This ratio deviates significantly and increasingly from 1 for <0.15 and A>1. Such moistures are quite common even at relatively humid mid-latitudes. From the fact that the v_d measured in the field can reach 0.1 cm/s (see Ehhalt and Rohrer, 2009), we can also conclude that A>1 in most natural soils and often reaches up to 8. That means for many field conditions we must expect the actual v_d to be smaller than that calculated from a one-layer model. At the same time, Fig. 8 indicates that there is no single correction formula to convert v_d from the one-layer model to a v_d calculated with the two-layer model. We, therefore, recommend the use of eq. (13) for calculating v_d, although it means a slightly larger calculation effort. This should lead to a significant improvement in the calculated v_d especially in dry climates.

Fig. 8. 

In fact, a general application of process-based calculations of v_d – whether based on a one-layer or two-layer model – is hampered far more by the lack of data on k_s(Θ_w,T)Θ_a. That moisture dependence has been determined only for three soils: In addition to the two soils described here, sand and loess loam, k_s(Θ_w,T)Θ_a has been determined for boreal forest soil by Smith-Downey et al. (2006). This limits a global calculation of soil uptake with the help of eq. (13). On the other hand, many soils consist of a mixture of loam and sand as the major components. Thus, appropriate interpolations between eqs. (16) and (17) for k_s(Θ_w,T)Θ_a and between eqs. (18) and (19) for δ could be used to treat those cases until measurements of k_s(Θ_w,T)Θ_a in further soils become available.

6. Summary and conclusions

We derived an analytical expression for the deposition velocity of H₂ on soil, v_d, which explicitly describes the dependency of v_d on the average soil moisture, , and temperature, T, of the soil, which are the only dependencies that have been identified so far. The derivation is based on a two-layer soil model originally proposed by Yonemura et al. (2000b). It consists of a dry top layer with a depth, δ, free of H₂ removal and a moist deeper layer with a uniform H₂ removal rate constant, k_sΘ_a. With a few further assumptions – notably uniform moisture and therefore uniform diffusivity, D_S, in each of the layers – we solve the 1D vertical diffusion equation in each layer. By matching the H₂ fluxes at the boundaries, we obtain three equations which can be solved for v_d in terms of δ, D_S, k_sΘ_a. For a quantitative estimate of v_d an approximate expression for δ() is derived. As it turns out δ is not a unique function of , but influenced by other factors such as the immediate history of rainfall and evaporation. The derived v_d shows a relatively weak dependence on T, but varies strongly with . The dependence on exhibits a broad maximum of v_d around =0.12 for loess loam and =0.07 for sandy soil. v_d approaches 0 as approaches 0.38, the value of Θ_P. The decrease of v_d to low is steep and reaches 0 at a threshold ~ 0.02 or 0.01, respectively. A comparison of the calculated v_d() with field measurements is hampered by the lack of suitable data. For >0.05, where field data on v_d exist, the agreement is reasonable. A comparison with numerically calculated v_d() from a completely characterized artificial model soil gives good agreement and favours the v_d calculated from the two-layer model over that calculated from a one-layer model.