Predicting the impact of spatial heterogeneity on microbial redox dynamics and nutrient cycling in the subsurface

Abstract. The subsurface is a temporally dynamic and spatially heterogeneous compartment of the Earth's Critical Zone, and biogeochemical transformations taking place in this compartment are crucial for the cycling of nutrients. The impact of spatial heterogeneity on such microbially mediated nutrient cycling is not well known which imposes a severe challenge in the prediction of in situ biogeochemical transformation rates and further of nutrient loading contributed by the groundwater to the surface water bodies. Therefore, we undertake a numerical modelling approach to evaluate the sensitivity of groundwater microbial biomass distribution and nutrient cycling to spatial heterogeneity in different scenarios accounting for various residence times. The model results gave us an insight into domain characteristics with respect to presence of oxic niches in predominantly anaerobic zones and vice versa depending on the extent of spatial heterogeneity and the flow regime. The obtained results show that microbial abundance, distribution, and activity are sensitive to the applied flow regime and that the mobile (i.e., observable by groundwater sampling) fraction of microbial biomass is a varying, yet only a small, fraction of the total biomass in a domain. Furthermore, spatial heterogeneity resulted in anaerobic niches in the domain and shifts of microbial biomass between active and inactive states. The lack of consideration of spatial heterogeneity, thus, can result in inaccurate estimation of microbial activity. In most cases this leads to an overestimation of nutrient removal (up to twice the actual amount) along a flow path. We conclude that the governing factors for evaluating this are the residence time of solutes and the Damköhler number (Da) of the biogeochemical reactions in the domain. We propose a relationship to scale the impact of spatial heterogeneity on nutrient removal governed by the log10Da. This relationship may be applied in upscaled descriptions of microbially mediated nutrient cycling dynamics in the subsurface thereby resulting in more accurate predictions of e.g., carbon and nitrogen cycling in groundwater over long periods at the catchment scale.


The reviewer is correct in assuming that it is arithmetic average. In addition, we would like to clarify that we ensured that the average water flux in all the domains in a particular flow regime was the same by adjusting the hydraulic gradient between the inlet and outlet boundaries (L167-169). We apologise for the confusion. To calculate the Peclet number, the domain length was used (0.5m) to derive the Peclet number for the flow processes occurring at the observation scale. To avoid confusion, we propose to delete this entry from Table 1.
Line 195: What is the correlation length used in the simulations? I couldn't find it in the tables.
We apologise for the confusion. The correlation length referred to in this sentence is the autocorrelation spatial/length scale. Its value is 0.1 m (20% of the domain size considering that several studies have found correlation length to be 10-30% of the observation scale (Turcke and Kueper (1996) found typical correlation length to be varying from 0.16m and 0.23m in core sizes of 1.5m while Welhan and Reed (1997) observed correlation length to be 1.5 km in the total domain size of 15 km)), We will update this sentence to prevent further confusion: " … Each random field was characterized by the same mean value of conductivity (i.e., average conditions at the subject site (Jing et al., 2017)) and spatial autocorrelation length scale (0.1 m) in all realizations, in scaling with the size of the domain in line with previous studies (Turcke and Kuper, 1996;Welhan and Reed, 1997;Desbarats and Bachu, 1994)."  The outcomes may be sensitive to the assumption of a second-order stationary random field with horizontal anisotropy. Other types of heterogeneity (e.g., multipoint statistical models, geometric models, or depositional process models) could lead to different (and probably even more striking) conclusions. I don't view this as a flaw of the study, since this assumption is conventional, but the assumption and its potential implications should perhaps be discussed.
We agree with the reviewer that other models are available to describe spatial heterogeneity. We selected the variogram method as it condensed the representation of the spatially heterogeneous fields in a limited parameter set (de Marsily et al., 2005). Thus, it was able to represent a wide variety of heterogeneous fields regardless of geological or depositional processes involved that resulted in the creation of such spatially heterogeneous domains in the first place, or without incorporating large datasets (as required in multipoint approaches). In our work, we overcame any bias induced by the variogram approach by using the breakthrough time to indicate the extent of spatial heterogeneity, independent of how heterogeneity is described. This led us to discuss results on impact on nutrient cycling and biomass in terms of reduction in breakthrough time given the same average flow conditions. We propose that the same approach could be followed regardless of the heterogeneity model used to generate such spatial random fields. We will add a note on this aspect in the Discussion section of the revised manuscript.
Lines 245-246: Clarification is needed here to elucidate what is meant by "fit of the model" in defining the AIC criterion. It isn't clear whether this is fit to actual field observations (there are some discussed in the paper but they are not described in detail) or fit to analytical solutions (e.g. Figure S8), or something else.
We understand the source of ambiguity, we have updated this section as follows: " To evaluate the key factors determining the impact of spatial heterogeneity on nutrient cycling, we undertook a series of multivariate statistical analyses of the simulation results using Linear Mixed Effect Modelling, progressively including variables in both fixed effects and random effects. We compared the Akaike Information Criterion (AIC) of each model to evaluate the fit of the model. AIC is an indicator of prediction error associated with a general linear model. It is an indicator of relative performance of a group of models; the model with the lowest AIC is concluded to be the one with least prediction error or best performance. With each iteration of the model, we selected the features most influencing the performance of the model and reducing the AIC of the predictions. " Line 253 (Equation 2): Is this a standard definition of Da? If so, please provide a citation. If not, please provide clarification as to how this equation represents the ratio of advective and reactive time scales.
We realize that this equation was not adequate to explain the derivation of Da in the scenarios. We will update the section with the following explanation on calculation of Da for the various scenarios. We also noticed there was an error in the equation which is rectified in the explanation below as well. " we used the Damköhler number (Da) to indicate the reaction regime for each reactive species. Da is defined as the ratio of the advective transport time scale and the reaction time scale as described in Eq. 2.
= , where, τreaction is the characteristic reaction time scale and τtransport is the characteristic transport time scale given by the breakthrough time of a conservative tracer in the domain. We adapted this definition and used Eq 3 below to calculate the apparent Da using values estimable in the field when > 5% (adapted from Pittroff et al., 2017).
where, 5 is the concentration of the chemical species at the first cross-section (y = y5) when ≤ 5%, , Thus, we were able to characterize reaction dominant system where Da > 1. We took the logarithm of Da to the base 10 (log10Da) to characterize the regime for each reactive species in each domain.

For a scalable relationship addressing impact of spatial heterogeneity on reactive species removal, we conduct a simple linear regression analysis of species removal vs. residence time (both in relative units to the homogeneous reference cases) for different log10Da ranges. "
Lines 265-270: These lines highlight a broader issue: What times were used for analysis and metrics of reactive processes? The flow field is clearly stated as being steady, but the concentration fields will be transient, and various times are mentioned in the manuscript. The breakthrough time is defined as the time at which Cout/Cin = 0.5 (for tracer), but what other times are considered for reactive species (they may not reach this value at the outlet)? This was confusing to me as a reader and should be clarified.
We agree that this aspect is confusing for the reader. We will rephrase the Data Analysis section to specify that (with the exception of the tracer tests) we are using the steady state concentration profiles for chemical and microbial species as well.  5). We used a wide range of reducing breakthrough times (normalized by that in base cases, from 10% to 90% of breakthrough time in the base case) to solve Equations 5 and 6 (subsequently plotted in Fig. S8.

Technical Corrections:
Lines 58-60: awkward wording, suggest rephrasing Thank you for the feedback. We will rephrase the sentence as follows: "Microbial communities play a key role in these biogeochemical cycles since they mediate nearly all the naturally occurring processes that contribute to these cycles." Lines 76-78: consider simplifying to "…representative of a system's chemical and biological species, and second…representative of a system's flow and transport pathways." Thank you for the feedback. We will rephrase the sentence as follows: "First, the reaction network should be representative of a systems' chemical and biological species, and second, the flow component of the model should be representative of a system's flow and transport pathways." Lines 102-103: suggest grouping citations at the end of the sentence

Noted and addressed.
Line 153: e.g. seems out of place, consider deleting

Noted and addressed.
Line 155: aerobic should be all lowercase

Noted and addressed.
Line 162: necromass is one word, not two

Noted and addressed.
Line 163: complete the second half of the sentence by describing how microbes become immobilized (biofilms etc.) Mathematically, the attachment of microbes depends on only the carrying capacity and concentration of immobile microbes (see sections A.3.3, L705 onwards). In a real system, we agree that immobilisation of microbes may be due to several reasons such as biofilms, interaction with the matrix etc. We rephrase the sentence as follows: "Furthermore, the reaction network accounts for microbial attachment, in case of hospitable conditions, and detachment due to inhospitable conditions or velocity of the water (see section A.3.3). The detached mobile bacteria are transported by the flowing water."

Addressed.
Line 270: this is "the" same (add "the" to the sentence)

Addressed.
Line 280: switch scale and spatial in the sentence Here, we intended to refer to the scale of the spatial scale sample. We will instead rephrase the sentence as follows for better readability: "We explore flux-averaged concentrations of mobile species and spatially averaged concentrations of immobile species in 1-D, along the predominant flow direction, and explore the 2-D concentration heat maps of the domain to compare the information lost when neglecting spatial heterogeneity at scales smaller than that of the sample." We agree, and we will remove the title in the revised manuscript.
Line 424: unbold sentence Addressed. This is a legacy error. We will remove the border in the revised figure.

Addressed
Line 482: Is there a difference between dormant and inactive? The term "dormant" is used sparingly within the manuscript, so I'd suggest sticking to inactive and defining that this pool includes dormant microbes to avoid confusion.
We agree that we should be consistent with the terminology. We will use "inactive" consistently throughout the revised manuscript.
Lines 487-489: clunky sentence, suggest rewording Rephrased in the revised manuscript: "The system may respond similarly to temporal fluctuations in groundwater velocities resulting from seasonal cycles as well." Figure 6: It's a bit unorthodox to present figures within the discussion section. Why didn't you introduce it within the results section and then reference it in the discussion?
We agree with the reviewer in that it is better to introduce Figure 6 in the Results section. We will rearrange the manuscript to accommodate this.
Line 618: Change to "The regression model links the..."