Spatial Control of Microbial Pesticide Degradation in Soil: A Model-Based Scenario Analysis

Microbial pesticide degraders are heterogeneously distributed in soil. Their spatial aggregation at the millimeter scale reduces the frequency of degrader–pesticide encounter and can introduce transport limitations to pesticide degradation. We simulated reactive pesticide transport in soil to investigate the fate of the widely used herbicide 4-chloro-2-methylphenoxyacetic acid (MCPA) in response to differently aggregated distributions of degrading microbes. Four scenarios were defined covering millimeter scale heterogeneity from homogeneous (pseudo-1D) to extremely heterogeneous degrader distributions and two precipitation scenarios with either continuous light rain or heavy rain events. Leaching from subsoils did not occur in any scenario. Within the topsoil, increasing spatial heterogeneity of microbial degraders reduced macroscopic degradation rates, increased MCPA leaching, and prolonged the persistence of residual MCPA. In heterogeneous scenarios, pesticide degradation was limited by the spatial separation of degrader and pesticide, which was quantified by the spatial covariance between MCPA and degraders. Heavy rain events temporarily lifted these transport constraints in heterogeneous scenarios and increased degradation rates. Our results indicate that the mild millimeter scale spatial heterogeneity of degraders typical for arable topsoil will have negligible consequences for the fate of MCPA, but strong clustering of degraders can delay pesticide degradation.


b. Reactive Transport
The advection-diffusion-reaction equation (ADR) was defined using COMSOL®'s Transport of Diluted Species in Porous Media module as given by Eq. 1 in the main text where [kg/m 3 ] is the soil bulk density and the relation between solution phase pesticide concentration [ mol C/m 3 ] and the sorbed phase concentration [ mol C/kg] was defined with the Freundlich sorption isotherm given by Eq. 2 in the main text as C S = K F (C L ) n F , with the Freundlich coefficient [( mol C/kg)( molC/m 3 ) − ] and exponent [1]. Spatial heterogeneity in the hydraulic parameters was not explicitly considered in our modeling framework and parameters were kept uniform in each depth layer.  (Vanderborght and Vereecken, 2007). To test the sensitivity of the model to the choice of , its value was varied to 0.01, 0.05 and 0.1 m and times needed for 50% degradation in the different heterogeneity scenarios were compared (supplementary Figure S4). We further assumed a ratio of = 3 but also tested ratios of = 10. A uniform value for both and was assigned to the entire simulation domain.
The soil diffusion coefficient [m 2 /s] accounts for decreased diffusivity under unsaturated conditions. This was accounted for in COMSOL® using the Millington and Quirk (1961) with the molecular diffusion coefficient of MCPA [m 2 /s] and the porosity ≔ .
The pesticide degradation rate [ mol C/m 3 /s] was given by Eq. 3 in the main text and was a function of microbial degrader concentration [ mol C/kg]. Microbes were considered to be in steady-state and non-mobile ( = 0). Degrader distributions were assigned via the Domain ODEs and ADEs module in COMSOL®. To allow for mass-balance closure a "degraded-C" pool ( 2 [ mol/kg]; i.e., assuming 100% of the degraded MCPA evolves as CO2) was considered as Mass balance in the entire soil column domain Ω (with the virtual thickness , [m]) was thus closed for MCPA as where ( = 0) [ mol C/kg] is the total initial MCPA concentration.
Initial and boundary conditions were assigned as described in the main text.

Scale transition theory
The second order accurate approximation of the reaction rate for Monod-type kinetics is given by Eq. 7 in the main text (Chakrawal et al., 2020) as where the partial derivatives are evaluated at the mean values of the variables ( ̅̅̅ and ̅ ) and ∑HOT represents the cumulative contribution from higher order terms. Note that the additional variance term 1 2 2 2 var( ) = 0 for Monod-type kinetics because the rate is linear with respect to (Chakrawal et al., 2020). The reaction rate ( , ) in this case is given by Eq. 3 in the main text as ( , ) = + . Consequently, the mean field approximation (MFA), the first right hand side term in Eq. S9 (Eq. 7 in the main text) is given as which is identical with the formulation in Chakrawal et al. (2020) except for needed here as a unit conversion factor, which carries over to the partial derivatives where ̅ and ̅̅̅ are the spatial averages ( ̅ ( ) = ∫ ∫ ( , , )d d ⋅ (∫ ∫ d d ) −1 ) of degrader and solution phase pesticide concentrations, respectively. Eq. S12 was used to assess the contribution of the covariance term to the difference between homogeneous and heterogeneous scenarios. As explained in the main text, the deviation from the mean field approximation (MFA); i.e., ̅ − MFA, was compared between homogeneous and heterogeneous scenarios.

Measurements of soil water retention curves, hydraulic conductivities and bulk densities
We sampled undisturbed 100 cm 3 soil cores from three soil layers (0-30, 30-60 and 60-90 cm; 20 replicates) by pushing stainless steel cylinders (5.6 cm diameter and 4 cm height) horizontally into the soil. The cylinders were then carefully excavated to avoid any loss of soil material on top and bottom surfaces. The soil cores were used to estimate soil water retention curves (5 replicates per depth) and saturated hydraulic conductivity (10 replicates) following standard procedures (DIN 19683-9; DIN EN ISO 11274). Bulk density was estimated by weighting oven-dried (105°C) soil cores (5 replicates) following the determination of soil water retention curves.

Estimation of soil hydraulic parameters
Measured soil water retention curves (SWRC) relate the pressure [Pa] (commonly given as pF value, = − log 10 ( 100 Pa ) ) to the respective water content. SWRC's were obtained for three depths in the reference soil and were individually fitted using the Mualem-van Genuchten model given in Eq. S3 reformulated for where the pressure head ℎ was computed from with the fluid density = 1 [kg/m 3 ] and the gravitational constant = 9.81 [m/s 2 ] as ℎ = . The saturated water content was directly obtained from the measurements and the hydraulic soil parameters [1], [1/m] and [1] were fitted with a hybrid global-local optimization algorithm with MATLAB®. First, MATLAB®'s particleswarm algorithm was used for global optimization and the found solution was then used with MATLAB®'s local fmincon opimization algorithm. particleswarm was run with a function S6 tolerance for termination of 1e-9 and otherwise default settings. fmincon was adapted to use the medium-scale sqp algorithm with optimality tolerance set to 1e-10 and 10,000 maximum function evaluations. Parameter bounds were 0-10 5 , 1-10 and 0-for , and , respectively. All measurement points, including replicates, were fitted simultaneously and the sum of squared residuals was minimized. Best fits to the experimentally measured SWRC are shown in Figure S1. Best fitting parameters are given in Table 1 in the main text.

Fitting of the MCPA degradation model
A process-based model that considers non-linear pesticide sorption and biodegradation was fitted to MCPA degradation data from recent microcosm experiments using the reference soil (Wirsching et al., 2020). The model accounted for Freundlich equilibrium sorption and Monod-type degradation kinetics by the following ordinary differential equation of the solution phase MCPA concentration: Diverging from the simulations with transient water content presented in the main text, experiments were run at a constant volumetric water content = 0.25 and soil bulk density , = 1200 kg/m 3 (Wirsching et al., 2020). Thus, equilibrium sorption was considered by the retardation factor = 1 + ( ) −1 , .
Wirsching et al. (2020) measured MCPA degradation by the indigenous microbial population in soil samples of the reference soil (an arable Luvisol; SM3, CAMPOS; for further details on soil properties and soil sampling see Wirsching et al. (2020)). In brief, tfdA gene abundance was measured as a proxy of the microbial degradation potential and degradation was assessed at different levels of initial MCPA concentration (0, 0.03, 0.05, 0.1, 0.5, 1, 5, 20 mg/kg). Since tfdA genes only slightly increased in response to highest amendment of the experiment, we used the average tfdA abundance (1.11 ⋅ 10 8 genes/kg) in the control treatment as a proxy for the MCPA degradation potential in the topsoil ( [µmol C / kg] = 1.11 ⋅ 10 8 genes/kg ⋅ / ). Timeseries of total residual pesticide concentration in the two batch experiments with the highest initial MCPA concentrations (5 and 20 mg/kg) were fitted with the same set of parameters. The differential equation in Eq. S14 was integrated with MATLAB®'s ode45 with relative and absolute tolerance set to 10 −10 and 10 −9 , respectively. As before for fitting SWRCs, all individual measurement points, including replicates, were fitted simultaneously. To better represent degradation dynamics at low concentration levels, all values were log10-transformed before computing sums of squared residuals. and where constrained to values between 10 −6 -10 3 and 10 −3 -10 9 , respectively. Fitting log10-transformed data was motivated by the better representation of degradation dynamics at low residual concentrations but caused initial degradation to be overestimated (compare supplementary Figure S2: original data (grey circles) vs. fit (red line)). This mismatch with experimental observations, however, does not impede inter-scenario comparisons.  (Wirsching et al., 2020). Data was fitted to log10-transformed residual concentration data (C,D).The black asterix in (A) and (C) mark data points below the limit of quantification (LOQ = 13 µg/kg; Wirsching et al., 2020)

Model exploration for different dispersivity values and ratios
Soil dispersivity values (Eq. S5) are highly uncertain (Vanderborght & Vereecken, 2007). We thus tested the sensitivity of half-lives (DT50, time when 50% of the initially applied pesticide dissipated from the entire simulation domain) to changes in from 0.01 to 0.1 m (Vanderborght & Vereecken, 2007) at ratios of 3 and 10 (default values were: = 0.03 m and / = 3, see Table 1 in the main text).
DT50 values in HOM and LOW hardly changed over the explored dispersivity range (supplementary Figure S5). Due to the increased mixing of substrate in EXTR (and to a lesser extend in HIGH), DT50 decreased with increasing dispersivities (both, and ). Intra-scenario variation of DT50 with the default dispersivity settings, however, largely exceeded the variation observed for the range of tested dispersivities in a single distribution realization (supplementary Figure S5 B, D). Figure S6 shows the same data as Figure 2 in the main text but with logarithmic y-axes.

Dimensionless scale transition approach
Wilson and Gerber (2021) recently suggested a dimensionless reformulation of the scale transition approach by Chakrawal et al. (2020). In brevity, Wilson and Gerber (2021, Eq. 17, p. 5672) suggest to express Eq. S9 (Eq. 7 in the main text) truncated to second order accuracy (i.e., neglecting ∑ ) in a dimensionless form which, applied to the Monod-type kinetics we assumed corresponds to (1 + ) 2 2 ) . (S17) These formulations allow for analytical examination of how ( ), respectively ( ) impact the systems deviation from the mean field approximation (i.e., the scale transition correction) for specific combinations of , , , and 2 values. In our simulations all these variables (with the exception of ( )) vary through time. How ( ) and ( ) influence the scale transition correction thus changes throughout the course of the simulation, and between heterogeneity scenarios. The analysis presented in the main text, based on the covariance term, accounts for these interactions as the covariance between degrader and substrate distribution can be expressed as ( , ) = , ( ) ( ) ̅̅̅̅ , (S18) similar to Eq. 9 in Wilson and Gerber (2021, p. 5671). Despite the multiple sources of variability, the dimensionless depiction of simulation outcomes illustrates that the relative deviation from the MFA (the scale transition correction) does not diminish as substrate is used up over time, but that a stable deviation is reached in the order HOM<LOW<HIGH<EXTR ( Figure S8).