Parameter sensitivity of SWAP–PEARL models for pesticide leaching in macroporous soils

Pesticide transport simulation by SWAP–PEARL (Soil–Water–Atmosphere–Plant and Pesticide Emission Assessment at Regional and Local scales) models can help to predict pesticide leaching at regional scales. For reasons of economic and time efficiency, measurement efforts should be prioritized towards critical parameters. The objective of this research is to perform a Morris screening and Sobol–Jansen sensitivity analysis to SWAP–PEARL models, using a reasonable worst‐case scenario. Three pesticide compounds were analyzed:, bentazon (zero sorption), imidacloprid (moderately sorbed), and compound I (highly sorbed). Initial macropore and pesticide parameter values were varied by ±20% to generate parameter ranges. The outputs analyzed were the concentration in drainage water, the average concentration in groundwater between 1 and 2 m, and the concentration in the soil system at 100‐cm depth. Influential parameters found through the Morris method were analyzed using the Sobol–Jansen method. The results for bentazon indicate that the degradation half‐life (DT50), the bottom depth of the internal catchment (zic), and the proportion of the internal catchment at the soil surface (pic_0) are critical parameters in all the outputs analyzed. For imidacloprid and compound I, the most relevant parameters for drainage output are the Freundlich sorption exponent (Fexp) and zic; for groundwater, the relevant parameters are Fexp, the bottom depth of static macropores (zst), and pic_0; and for soil concentrations at 100‐cm depth, the relevant parameters are Fexp, zic, and pic_0. The Morris and Sobol–Jansen methods produce the same results for the first position in the ranking. Measurement efforts should be performed to update national soil databases, including critical pesticide and macropore parameters.


INTRODUCTION
The use of plant protection products in agriculture can result in leaching of these substances to groundwater and emissions into surface water via drainage. These environmental risks need to be assessed before these products can be used in agriculture. The procedure for authorization in the United States was developed by the USEPA (2020). In the European Union, the assessment must comply with EU Regulation 1107/2009. To assess the risk of leaching to groundwater and the risks to aquatic organisms, groundwater and surface water scenarios have been developed: FOCUS (2000FOCUS ( , 2001 and EC (2014). In the FOCUS groundwater scenarios, the MACRO model (Jarvis & Larsbo, 2012) is applied to calculate the transport of substances in macroporous soils. Other models that can describe the water fluxes in macroporous soils are the Soil-Water-Atmosphere-Plant model (SWAP), in which a dualpermeability concept has been implemented (Kroes et al., 2017), and HYDRUS (Šimůnek, van Genuchten, & Šejna, 2016). The SWAP model has been coupled to the pesticide fate model PEARL (Pesticide Emission Assessment at Regional and Local scales; van den Berg, Tiktak, Boesten, & van der Linden, 2016). The combination of these models is the kernel of FOCUSPEARL. This model also facilitates the utilization of the European FOCUS groundwater scenarios.
Various parameters that describe the characteristics of the macropore system have been added to SWAP-PEARL (Tiktak, Hendriks, Boesten, & van der Linden, 2012;Tiktak, Hendriks, & Boesten, 2012). For example, considering a single soil horizon in a macroporous soil, at least 15 and 7 parameters are required in the SWAP and PEARL models to describe the hydrology and behavior of pesticides, respectively. Those parameters should be measured, preferably onsite, to assess the behavior of the pesticide at the field scale, which may be time consuming and expensive. Therefore, the measurement effort should focus on the most critical parameters involved in simulating pesticide concentrations in the soil, groundwater, and drainage systems. A ranking of parameter importance may be constructed through a sensitivity analysis.
A sensitivity analysis is defined as the study of how uncertainty in the model output can be apportioned to different sources of uncertainty in the model input (Saltelli et al., 2010). Saltelli (2002) defines factor prioritization as the use of sensitivity analysis to discover the most relevant parameters regarding a model output. He also described factor fixing as the use of sensitivity analyses to find nonimportant parameters concerning a model output. Global sensitivity analyses are commonly used for nonlinear and nonadditive models (Gan et al., 2014;Pianosi et al., 2016;Saltelli et al., 2008;Song et al., 2015). Global sensitivity analyses can be qualitative, such as the Morris screening method (Morris, 1991), or quantitative, such as the Sobol-Jansen method (Saltelli et al., 2010). Both the Morris and Sobol-Jansen methods have been combined to analyze complex hydrological models (Cuntz et al., 2015).

Core Ideas
• Two global sensitivity analyses were performed on pesticide leaching in macroporous soils. • The Morris and Sobol-Jansen methods yield comparable parameter importance rankings. • Sorption and degradation are key factors when assessing leaching in macroporous soils. • Depths of static macropores and internal catchments are also critical soil parameters. • The addition of macropore parameters in databases of soil spatial data is recommended.
A sensitivity analysis for SWAP-PEARL has been performed for homogeneous soils (Boesten, 1991), but there is little data on the sensitivity of hydrological and substance fluxes, and concentrations in macroporous soils, to changes in input parameters. For the current sensitivity analysis, Andelst soil was selected (Scorza Júnior, Smelt, Boesten, Hendriks, & van der Zee, 2004;Tiktak, Hendriks, Boesten, & van der Linden, 2012;Tiktak, Hendriks, & Boesten, 2012). This soil contains a macropore system with dead-end macropores, structural and shrinkage macropores, tile drains, and close interaction between the subsoil and the shallow groundwater. In Andelst soil, fast transport of bentazon and imidacloprid to drainage and groundwater systems was detected due to macropore flow (Scorza Júnior et al., 2004). The strength of macropore flow in the Andelst soil is also demonstrated in Scorza Júnior and Boesten (2005), where the chromatographic version of SWAP-PEARL could not simulate the pesticide leaching for both compounds. Therefore, macropore flow is the dominant process in the Andelst soil, making it a suitable study case for sensitivity analysis.
The objective of this research was to perform a Morris screening and Sobol-Jansen sensitivity analysis of the SWAP-PEARL models using Andelst soil as a reasonable worst-case scenario. The outcome of this research can be used as a guide for prioritizing macropore and pesticide parameters to be included in national soil databases for the three contrasting pesticides selected.

Morris elementary effect screening method
A factor in sensitivity analysis is any model input that can produce variations in the model output (Saltelli et al., 2008). In this paper, the term factor is used interchangeably with the term parameter. The Morris elementary effect screening method is a global sensitivity one-at-atime method and constitutes a significant improvement of the gradient-based method (Morris, 1991;Pianosi et al., 2016;Song et al., 2015). The methodology computes elementary effects (EE i ) for the factor X i by generating trajectories (r) in the input factor space (Equations 1a and 1b): where Y is the model output, Δ is the change on the factor X i calculated by Equation 1b, k is the total number of factors, p is the number of levels that factor X i is allowed to vary in the parameter space, and ω is a scalar called "grid jump" commonly estimated as p/2 (Morris, 1991;Pujol, 2017). The generation of trajectories (r) can be seen in Saltelli et al. (2008). Cuntz et al. (2015) mentioned that the number of trajectories (r) should be at least equal to the number of factors. Both the mean (μ) and the standard deviation (σ) are computed over all trajectories calculated for the factor X i : The overall effect of X i in the output variation is μ i , whereas σ i represents the interaction between factors. The main limitation of Equation 2 is that positive and negative EE i realizations cancel each other out, incorrectly indicating that the factor is not essential (Type II error in Saltelli et al., 2008). Campolongo, Cariboni, and Saltelli (2007) proposed the use of Equation 4 to overcome that issue: where μ * is the modified overall effect introduced by Campolongo et al. (2007). The statistics calculated by the Morris elementary effect method (μ i , μ * , and σ i ) generate a ranking of factor importance (Song et al., 2015), but not their magnitudes. Therefore, the Morris method is qualitative.

Sobol-Jansen sensitivity analysis
The quantification of factor importance is commonly performed using variance-based methods (Song et al., 2015). Those methods are based on the law of total variances (see Yun, Lu, & Jiang, 2018). The main effect, S i , is commonly used for factor prioritization, whereas the total effect, S Ti , is used for factor fixing (Saltelli et al., 2008). Saltelli et al. (2010) proposed the estimation of S i and S Ti as follows: where V[f(A)] is the total variance of model outputs computed from matrix A, A and B are two matrices of dimension (n s , k), n s is the number of rows in each matrix A and B that is related to the number of samplings points, and k is the number of columns in each matrix which is equal to the total number of factors. f(B) j , [ ( ) ] , and f(A) j are the model output for parameters in row j of the matrices B, ( ) , and A, respectively, and ( ) is a matrix with all the columns from matrix A except the ith column, which comes from matrix B. Therefore, the matrix ( ) is generated for each factor. More details about the construction of matrices A and B can be found in Cuntz et al. (2015). Equation 6 is based on the equation used by Jansen (1999), and the sampling method for filling matrices A and B in Equations 5 and 6 is the Sobol procedure. Therefore, this sensitivity analysis method is customarily called the Sobol-Jansen method. An n s value >1,000 is recommended for the convergence of S i and S Ti in Equations 5 and 6 (Cuntz et al., 2015;Gan et al., 2014;Pianosi et al., 2016).

Field site
The macroporous soil in Andelst (the Netherlands, Gelderland Province) was chosen for performing the sensitivity analysis. The Andelst study has been described extensively elsewhere ( The top 3 m of soil in Andelst is clay, and below that, sandy layers are present. Cylindrical macropores were observed in the soil down to 100-cm depth, arising from root decay and earthworm activity. Therefore, a mixture of rectangular (shrinkage cracks) and cylindrical macropore shapes (biopores) was identified. Tile drains were installed in the field site at 80-to-90-cm depth with a spacing of 1,000 cm. A winter wheat (Triticum aestivum L.) crop was cultivated twice during the years 1997 and 1999. Two pesticides with contrasting properties, bentazon and imidacloprid, were applied. Meteorological data, soil temperature, groundwater levels, drain discharge, and pesticide concentrations in drainage, groundwater, and the soil system were measured. The upward and downward water flow between a second aquifer and groundwater was determined from piezometer measurements at different depths. A graphic of the Andelst site is depicted in Figure 1.

SWAP model description
SWAP version 4.0.8 was selected for simulating water flow in the soil-plant-atmosphere continuum (Kroes et al., 2017). SWAP includes a multidomain approach to represent macropores (Figure 1, macropore system). The internal catchment domain is a set of dead-end macropores that are not connected to each other. The main bypass domain is a set of interconnected macropores that penetrate deep into the soil domain generating rapid drainage. The multidomain approach implies that the total relative macroporosity (w f ) changes over the soil profile. The meaning of the total relative macroporosity in SWAP is the sum of the relative macroporosity of the internal catchment and the main bypass domains over depth. The change of w f over depth is calculated mainly with the following parameters: the total relative macroporosity at the soil surface (w f_0 , cm 3 cm −3 ), the proportion of internal catchment at the soil surface (p ic_0 ), the bottom depth of static macropores (z st , cm), a shape factor parameter for the w f curve (m), and the bottom depth of internal catchment (z ic , cm). The analytical equations for representing the multidomain approach in SWAP are shown in the supplemental material as Supplemental Equations S2-S4. SWAP also considers dynamic macropores as produced by shrinkage cracks. The dynamic macropore volume is added to the corresponding macropore domain, and it is included in the total relative macroporosity. In this sensitivity analysis, dynamic macropores are simulated, TA B L E 1 The depth of the soil horizons, the van Genuchten-Mualem matrix parameters (residual water content [θ r ], saturated water content [θ s ], inverse air-entry value [α], pore size distribution index [n], saturated hydraulic conductivity of the matrix [K s ], and pore connectivity [l]), the entry pressure head (h w ), the bulk density (D a ), and organic matter (OM) for the Andelst soil (Scorza Júnior et al., 2004;Tiktak, Hendriks, Boesten, & van der Linden, 2012;Tiktak, Hendriks, & Boesten, 2012) Depth Three parameters mainly define the lateral exchange of water between the matrix and macropores. These are the maximum and minimum polygon diameter (d polmax and d polmin , respectively; cm), an empirical parameter for modifying the Parlange (1975) sorptivity equation (S Parlange ), and a practical parameter for adjusting the Darcy equation (S Darcy ). S Parlange and S Darcy are used in SWAP to calculate lateral exchange of water under dry and wet conditions, respectively. The equations for describing these processes can be found in the supplemental material (Supplemental Equations S5-S7).
The top boundary condition allows precipitation, rainfall interception by plants, and runoff to be included. The bottom boundary condition models seepage flux between the second aquifer and the groundwater considering aquitard resistance ( Figure 1). A lateral boundary condition is included to simulate drainage systems.
The drainage flux is computed for macropores and matrix domains separately. The rapid drainage parameter (R drares ) is the drainage resistance in the main bypass macropore domain. In contrast, the drainage resistance (D rares ) is used for the matrix. The equations can be found as Supplemental Equations S8 and S9. The description of the SWAP model indicates that it is capable of realistically simulating the hydrology of the Andelst site.

SWAP parametrization in Andelst
The winter wheat crop was simulated using the simple crop option in SWAP. The initial water content of the soil profile was set at hydrostatic equilibrium with the initial groundwater level at 81-cm depth. The potential evaporation and transpiration were computed using the refer-ence Penman-Monteith evapotranspiration obtained from meteorological data. Rainfall intensities were measured at the site. The soil profile was simulated down to 320-cm depth ( Table 1). The residual water content (θ r , cm 3 cm −3 ), the saturated water content (θ s , cm 3 cm −3 ), the inverse of the air entry value (α, cm −1 ), the pore size distribution index (n), the saturated hydraulic conductivity of the matrix (K s , cm d −1 ), the pore connectivity (l), the entry pressure head (h w , cm), the bulk density (D a , g cm −3 ), and organic matter (OM, kg kg -1 ) were estimated in the laboratory and are presented in Table 1. An atmospheric top boundary condition with surface runoff was included in the SWAP simulations. At the bottom boundary, the flux is computed from the hydraulic head of the deep aquifer and the simulated groundwater level. A lateral boundary condition includes drainage flow, which was simulated using the hydraulic head difference between drainage and surface water level and a drainage resistance. Heat flow parameters were obtained using measurements of soil texture and organic matter. The macropore system includes static and dynamic macropores. Macropore parameters were obtained through field observations and pedotransfer functions (Scorza Júnior et al., 2004;Tiktak, Hendriks, Boesten, & van der Linden, 2012;Tiktak, Hendriks, & Boesten, 2012). The parameters for the clay shrinkage curve were obtained under laboratory conditions. The initial macropore parameters and drainage resistance set in SWAP can be found in Table 2. 3.3 PEARL model

PEARL model description
The PEARL model is used to simulate pesticide behavior in soils with a macropore system. Using the soil hydrology TA B L E 2 The macropore parameters (z ic , z st , w f , p ic , m, d polmin , d polmax , S Parlange , S Darcy , and R drares ) and drainage resistance (D rares ) estimated for the Andelst soil (initial). Minimum and maximum values were computed as a 20% variation of the initial value computed by SWAP, PEARL solves the mass conservation equations for pesticides in soil.
The lateral exchange of pesticides between the matrix and macropore domains is by convection only, and the SWAP macropore parameter d pol is involved in this process. The transformation rate of pesticides is calculated with the degradation half-life parameter (DT50, d, Supplemental Equation S11). The sorption of pesticides to the soil matrix and macropores (Supplemental Equations S12 and S13) is computed with the Freundlich equation, using the following parameters: the profile of the soil organic matter fraction over depth, the organic matter/water distribution coefficient for the equilibrium sites (K OM , m 3 kg −1 ), and the Freundlich exponent (F exp [-]). Adsorption in macropores is only considered in the main bypass domain, including an average mass fraction of organic matter over the depth fraction of the water-filled main bypass domain. The parameter used for adsorption in the main bypass domain was set fixed to 0.02 as in previous studies (Tiktak, Hendriks, & Boesten, 2012) and is therefore not included in the sensitivity analysis.
Pesticide uptake by plant roots depends on the transpiration stream concentration parameter (P uptk ). The dispersion coefficient is computed using the dispersion length parameter (d L , m). Diffusion is described by Fick's law, using the reference diffusion coefficient in water (D dif , m 2 d −1 ).
The pesticide discharge by drainage from the matrix is calculated from the volumetric flow rate of drain water by multiplying it with the concentration of pesticide in the liquid phase. Drainage by macropores is calculated for the main bypass domain by multiplying the volumetric flow rate of water with the pesticide concentration in the main bypass domain. Finally, the concentration in drainage water is considered to be a mix between matrix and rapid drainage, weighted by their respective volumetric fluxes of water. A more detailed description of the model is found in the supplemental material.

PEARL parametrization in Andelst
A research version of the PEARL model version 3.2.2 was used to simulate the pesticides, bentazon, imidacloprid, and compound I. The last pesticide was not applied in Andelst, but it is used in the FOCUS surface water scenarios. This substance has been included to cover a broader range of pesticide properties.
Bentazon is very mobile in the soil with low persistence and negligible sorption, imidacloprid is moderately sorbed and persistent in the soil, and compound I is strongly sorbing and very persistent in the soil. The parameters for bentazon and imidacloprid were obtained from the initial measurements reported at the Andelst soil. The parameters for the compound I were obtained from the FOCUS surface water report (FOCUS, 2001).
The initial pesticide concentrations in the soil were set to zero. The three compounds were applied on the same date, 7 Apr. 1998, to make the results easier to interpret. The application rates of bentazon, imidacloprid, and compound I at the soil surface were 1.4, 0.55, and 0.1 kg ha −1 , respectively. Tillage was applied on 8 Dec. 1998 to a depth of 27.5 cm. The parameters for simulating the pesticides in the PEARL model are listed in Table 3.
In Tables 2 and 3, only parameters investigated in this sensitivity analysis are listed.

Morris elementary effect screening method included
The Morris elementary effect screening method was used to reduce the number of input factors for the Sobol-Jansen method. The SWAP macropore parameters incorporated in the sensitivity analysis are listed in Table 2. The soil hydraulic parameters of the matrix (Table 1) were not included.
The PEARL parameters incorporated in the sensitivity analysis (Table 3) for bentazon include d L , DT50, P uptk , and D dif . The parameters for imidacloprid and compound I are the same as for bentazone plus F exp and K OM . The dispersion length d L was taken to be constant over the soil profile for all the compounds.
The initial values of the parameters included in the Morris screening method were varied by ±20% to generate lower (minimum) and upper (maximum) parameter ranges (Tables 2 and 3). The percentage variation used to generate the parameter range was selected by trial and error, under the condition that all the model simulations with the SWAP-PEARL model converged in both sensitivity analyses. The Morris elementary effect parameters for bentazon were r = 100, k = 15, and Δ = 0.6, where p = 6 and ω = 3 (Equations 1-4). The Morris elementary effect factors for imidacloprid and compound I are the same as for bentazon, except for k, the values being 17 and 16, respectively. The Morris sampler generator and sensitivity index com-putation were performed with the R package "Sensitivity" (Pujol, 2017).
The outputs of the PEARL model analyzed include the concentration of pesticide in the drainage water from matrix and macropores ("drainage"), the average concentration of pesticide in the groundwater between a depth of 1 and 2 m ("groundwater"), and the concentration of pesticide in the soil system at 100-cm depth ("soil at 100-cm depth"). The PEARL outputs covered a period from 1 Jan. 1998 until 26 Apr. 1999. The sensitivity indices σ (Equation 3) and μ* (Equation 4) were computed for each day.
The parameter importance was evaluated following the classification system of Lammoglia et al. (2017). They defined three ranges of importance for Morris sensitivity analysis; highly influential if μ* > 0.5 μ * max , influential if 0.1 μ * max < μ* < 0.5 μ * max and low impact if μ* < 0.1 μ * max . Even though daily values of the sensitivity indices were obtained, we analyzed the parameter importance from the average of the sensitivity indices over time and at the arrival and peak of the pesticide compound.

The Sobol-Jansen method
The Sobol-Jansen sensitivity analysis method was applied to the highly influential and influential parameters for each pesticide and output from the Morris method.

F I G U R E 2 Average Morris sensitivity indices σ and μ* and Morris sensitivity indices at arrival and peak for drainage concentration.
The names of the three pesticide compounds are shown on the right-hand side of the graph. The variables are degradation half-life (DT50), Freundlich sorption exponent (F exp ), equilibrium sorption coefficient on organic matter (K OM ), dispersion length (d L ), plant uptake (P uptk ), the proportion of the internal catchment at the soil surface (p ic_0 ), depth of internal catchment (z ic ), relative macroporosity at the soil surface (w f_0 ), minimum polygon diameter (d polmin ), maximum polygon diameter (d polmax ), Darcy factor related to lateral water exchange (S Darcy ), depth of static macropores (z st ), shape factor (m), reference rapid drainage resistance (R drares ), and drainage resistance (D rares ) The number of sampling points for each factor in the Sobol-Jansen method was set as n s = 2,000. The Sobol sequence, including scrambling by Owen and Faure-Tezuka, was generated using the R package "Randtoolbox" (Petr, 2018). Scrambling further homogenizes the distribution of points in parameter space for quasi-random sequences in high dimensions (Dimov, Georgieva, Ostromsky, & Zlatev, 2013), and it has been demonstrated that this improves the Sobol method compared to unscrambled sequences (Chi, Beerli, Evans, & Mascagni, 2005). The sensitivity indices, main and total effect, and their 95% confidence intervals obtained through bootstrapping (n boot = 1,000) were computed using the R package "Sensitivity" (Pujol, 2017).
The parameter range used in the Sobol-Jansen method is the same as that applied in the Morris method (Tables 2  and 3). In the Sobol-Jansen method, the same outputs as with the Morris method were analyzed. The main and total effect was computed for each day of the time series. The days with zero pesticide concentration were discarded because the sensitivity indexes are zero.
A convergence criterion for the Sobol-Jansen method is that the sum of the main effect should not be higher than 1.0, and the sum of the total effect should not be less than 1.0. We discarded daily values that did not converge, applying a threshold of 1.1 for the main effect and 0.9 for the total effect.

Morris elementary effect screening method
Only drainage outcomes are presented in Figure 2. The average concentration in the groundwater between depths of 1 and 2 m and the concentration in the soil system at 100-cm depth are depicted in Supplemental Figures S1 and S2, respectively.
In red are the parameters classified as highly influential and influential for the three outputs analyzed (Figure 2, Supplemental Figures S1 and S2). They were combined for each compound to be used in the Sobol-Jansen method. For bentazon, the parameters selected were DT50, P uptk , p ic , z ic , d L , d polmax , w f , S Darcy , z st , and m. For imidacloprid, the chosen parameters were F exp , z ic , p ic , K OM , DT50, m, d polmax , d polmin , w f , S Darcy , z st , D rares , R drares , and d L . For compound I, the parameters selected were F exp , z ic , K OM , p ic , m, w f , d polmax , d polmin , and z st .

Sobol-Jansen method
The daily outcomes of the main effect are depicted in Figure 3. It can be observed that, generally, they are not constant over time.

F I G U R E 3
Daily outcomes of Sobol-Jansen main effect for the concentration of pesticide in the drainage flux from matrix and macropores (drainage), the average concentration of pesticide in the groundwater between depths 1 and 2 m (groundwater), and the concentration of pesticide in the soil system at 100-cm depth (soil at 100-cm depth) outputs. The dashed black line depicts the concentration of each compound over time, as simulated by SWAP-PEARL (secondary y axis). The variables are degradation half-life (DT50), Freundlich sorption exponent (F exp ), equilibrium sorption coefficient on organic matter (K OM ), dispersion length (d L ), plant uptake (P uptk ), the proportion of the internal catchment at the soil surface (p ic_0 ), depth of internal catchment (z ic ), relative macroporosity at the soil surface (w f_0 ), minimum polygon diameter (d polmin ), maximum polygon diameter (d polmax ), Darcy factor related to lateral water exchange (S Darcy ), depth of static macropores (z st ), shape factor (m), reference rapid drainage resistance (R drares ), and drainage resistance (D rares ) In general, the interaction between parameters was low. This outcome is observed through a comparison of the main and total effect (Figure 4). If they are close, interaction is low; otherwise, it is high. We report high interaction for the concentration of imidacloprid and compound I in the soil system at 100-cm depth (Figure 4).

Main effect for the pesticide concentration in drainage water
The main parameters that explain the concentration of bentazon in drainage water are z ic and DT50, whereas for the peak concentration in drainage, they are DT50 and p ic_0 . For imidacloprid and compound I, the main parameters are F exp and z ic (see Figure 3).

Main effect for the pesticide concentration in groundwater
For the average concentration of bentazon in the groundwater between depths of 1 and 2 m, the macropore param-eters z ic and p ic_0 are the most important, with some relevance of DT50 at the end of the period. For imidacloprid, the main parameters are F exp , z st , and p ic_0 . We observed that just before the peaks of concentration, the importance of z st is very close to that of F exp (Figure 3). For compound I, the situation is similar to imidacloprid; however, z st is more critical than F exp just before the concentration peaks.

Main effect for pesticide concentration in the soil system at 100-cm depth
For the concentration of bentazon in the soil system at 100-cm depth, the main parameters are z ic , p ic_0, and DT50. For imidacloprid, the most relevant parameters were z ic , p ic_0 , and F exp . It is observed in Figure 3 that z ic was the most relevant parameter just before the rate of imidacloprid concentration increase. For compound I, the results were the same as those obtained for imidacloprid.

F I G U R E 4
Average outcomes of the Sobol-Jansen main and total effect for the concentration of pesticide in the drainage water from matrix and macropores (drainage), the average concentration of pesticide in the groundwater between depths of 1 and 2 m (groundwater), and the concentration of pesticide in the soil system at 100-cm depth (soil at 100-cm depth) for the three pesticide compounds (right). The variables are degradation half-life (DT50), Freundlich sorption exponent (F exp ), equilibrium sorption coefficient on organic matter (K OM ), dispersion length (d L ), plant uptake (P uptk ), the proportion of the internal catchment at the soil surface (p ic_0 ), depth of internal catchment (z ic ), relative macroporosity at the soil surface (w f_0 ), minimum polygon diameter (d polmin ), maximum polygon diameter (d polmax ), Darcy factor related to lateral water exchange (S Darcy ), depth of static macropores (z st ), shape factor (m), reference rapid drainage resistance (R drares ), and drainage resistance (D rares )

Convergence of Sobol-Jansen method
The convergence of the Sobol-Jansen method was successful; we eliminated for the concentration of bentazon in drainage output 8 out of 213 daily points, yielding a discard percentage of 3.75%. For groundwater and soil at 100-cm depth outcomes, the discard percentage was 0%. For imidacloprid, the discard fractions for drainage, groundwater, and soil at 100-cm depth outcomes were 7.94, 0.8, and 1.86%, respectively. For compound I, the discard fractions for drainage, groundwater, and soil at 100-cm depth outcomes were 9.43, 1.06, and 1.05%.

Comparison between Morris screening and Sobol-Jansen methods
The comparison of the sensitivity analysis methods is performed by utilizing the average information of μ* for the Morris method and the average main effect for the Sobol-Jansen method (see Figures 2 and 4, Supplemental Figures S1 and S2).
For the parameter in the first position, the match between the Morris and Sobol-Jansen methods is perfect for all pesticide compounds and outputs. For the second position, the match is generally good except for the concentration of bentazon in the soil system at 100-cm depth, where the second and third position shifts. A similar situation occurs for the compound I in groundwater.

Critical parameters for pesticide leaching
The quantification of parameter importance by the Sobol-Jansen method using the main effect in the time series (Figure 3) indicates for bentazon that DT50, z ic , and p ic_0 are essential parameters in all the outputs analyzed. For imidacloprid and compound I, the most relevant parameters for the drainage output are F exp and z ic ; for groundwater, the relevant parameters are F exp , z st , and p ic_0, and for soil at 100-cm depth, the relevant parameters are F exp , z ic , and p ic_0 .
Previous research has found that the concentration of pesticides in the groundwater was more sensitive to F exp than to K OM (Tiktak, Swartjes, Sanders, & Janssen, 1994). Heuvelink, Burgers, Tiktak, and van den Berg (2010) found that DT50 was critical for the average pesticide leaching concentration at 1-m depth computed with the chromatographic version of the PEARL model. Similarly, Scorza Júnior and da Silva (2011) analyzed the cumulative mass of leached pesticide at 1-m depth for a very mobile compound and a moderately sorbed compound. They found DT50 and F exp to be critical for the three contrasting soil types from the Dourados river watershed, including heavy clay soil. van den Berg et al. (2012) concluded that DT50 and F exp were the most relevant parameters for the predicted environmental concentration at 1-m depth in soils without macropore systems using the GeoPEARL model. Dubus and Brown (2002) and Dubus, Brown, and Beulke (2003) suggested that in the presence of macropore flow, compound properties such as DT50 and F exp still exert a significant influence on leaching when using the MACRO model (Jarvis & Larsbo, 2012).
The results of this study, which include a macropore system simulated by the SWAP-PEARL models, agree with the studies mentioned above regarding pesticide-related parameters (DT50 and F exp ). However, this research shows that macropore parameters related to geometry are also relevant for pesticide leaching, especially at earlier times ( Figure 3). Therefore, macropore parameters z ic , p ic_0 , and z st are very important for the initial breakthrough of pesticides. For bentazon, the initial breakthrough is mainly explained by macropore parameters in all the outputs (Figure 3). For the moderately sorbed compound imidacloprid and the highly sorbed compound I, F exp exerts a large influence over the whole period (Figure 3). At early times, we observed higher importance of macropore parameters, but that vanished over time.
The macropore parameters z ic , p ic_0 , and z st are mathematically related to the variation of the relative macroporosity over depth (Supplemental Equations S2-S4). Therefore, advances in field or laboratory determination of the variation of the relative macroporosity over depth, along with DT50 and F exp , should be a priority for national scale databases.
The measurement of DT50 and F exp is mainly performed under laboratory conditions (Scorza Júnior et al., 2004), including some adjustments for field conditions (Tiktak, Hendriks, & Boesten, 2012). The geometrical macropore parameters are more difficult to obtain. van Schaik, Hendriks, and van Dam (2010) obtained some macropore geo-metric parameters by inverse estimation, including w f_0 and p ic_0 , whereas direct parameter estimation was performed for z st , d polmax , and d polmin , using dye staining. Urbina, van Dam, van den Berg, Ritsema, and Tang (2020) developed a methodology to obtain the relative macroporosity and d pol with disk infiltrometers. This methodology has the potential to be used for commonly encountered macropore shapes, and it can be applied at different cross-sectional areas over depth. Some pedotransfer functions are shown in Tiktak, Hendriks, Boesten, and van der Linden (2012) to apply SWAP-PEARL at a national scale. Although the parametrization of SWAP-PEARL remains a topic of discussion and development, previous studies can be used for a direct parametrization of the model under field conditions, including the most relevant parameters detected in this study.

Framing errors
Framing errors mentioned in Saltelli et al. (2008) refer to choosing the wrong parameter ranges or not including parameters that may be relevant in the sensitivity analysis. The modification of parameter values by 20% might seem unrealistic and a source of framing error (Tables 2 and 3). However, this parameter range was selected by trial and error with the goal of not producing gaps in the sensitivity analysis. Exploring the full parameter space is challenging in mechanistic models because unrealistic parameter combinations can be generated, resulting in nonconvergence of the numerical solution. Such gaps might be a problem when using the Sobol-Jansen method, given the uniformity of the quasi-random sequence of the Sobol sampling method (Saltelli et al., 2010). The exclusion of soil matrix parameters from the SWAP model might be another source of framing error because they can interact with macropore or pesticide parameters, modifying the final parameter importance. We nevertheless left them out to increase the computational efficiency of the SWAP model. Additionally, previous studies considering only matrix systems with the SWAP-PEARL models suggest minor importance of matrix parameters regarding pesticide parameters in pesticide outputs (Heuvelink et al., 2010;Scorza Júnior & da Silva, 2011).

Main effect
The main effect over time was used to determine parameter importance. The interpretation of the main effect can be difficult when interactions between parameters are strong (Pianosi et al., 2016). Our results indicate that interactions were not critical for both pesticides and all outputs, except for the concentration in the soil system at 100 cm for imidacloprid and compound I, for which interaction was considerable ( Figure 4). The interactions between parameters are reduced by decreasing the number of parameters involved in the Sobol-Jansen method (Song et al., 2015), which we try to accomplish by performing first the Morris screening method discarding noninfluential parameters. One central assumption in the Sobol-Jansen sensitivity analysis method is that the parameters are independent. Therefore, correlations between parameters might lead to erroneous interpretations of the main effect (Song et al., 2015). Kucherenko, Tarantola, and Annoni (2012) generalized the Sobol sensitivity analysis for correlated factors. They found that depending on the level of correlation, the main effect can be higher than the total effect, S Ti . Pianosi et al. (2016) indicated that the value of S Ti tends toward zero for correlations close to unity between parameters. Taking the above studies into account, a strong correlation among the parameters included in the Sobol-Jansen method was not observed. Therefore, it can be stated that the parameter importance computed by the main effect is adequate and free of strong interactions and correlations among parameters for the majority of outputs under analysis.

Comparison between Morris and Sobol-Jansen methods
The Morris and Sobol-Jansen methods agree on the most critical parameter for both pesticides, and all the outputs analyzed. They are different in the second position for two outputs. This outcome is an important verification step because two different global sensitivity analyses provide identical parameter importance for the first position and generally similar importance for the second position. Campolongo et al. (2007) found that the statistics μ * (Equation 4) obtained from Morris method is a good proxy for the total effect, S Ti (Equation 6). In this research, the general low interaction between parameters results in a good match between the main and total effect for the first positions. This might explain the perfect match between the Morris and Sobol-Jansen methods.
The fact that both rankings are generally comparable is essential because the Morris screening method requires fewer iterations than the Sobol-Jansen method. As an example, for imidacloprid, the number of iterations for the Morris method was 1,800, whereas the Sobol-Jansen method was 32,000. This outcome should be considered with care given the limited scenarios presented in this research. Therefore, further research should be performed to study the similarities of the rankings determined with both methods.

CONCLUSIONS
Morris elementary effect screening and Sobol-Jansen sensitivity analysis methods were applied to SWAP-PEARL models. Andelst soil was chosen as the input soil for both methods because it involves complex hydrological processes and macropore flow. Three pesticides were analyzed: bentazon (zero sorption), imidacloprid (moderately sorbed), and compound I (strongly sorbed). The outputs from the PEARL model, with which the sensitivity analyses were performed, were the flux-averaged concentration of pesticide in the drain water, the average concentration of pesticide in the groundwater between depths of 1 and 2 m, and the concentration of pesticide in the soil system at 100-cm depth. For the nonsorbing compound bentazon, the degradation half-life (DT50), the bottom depth of internal catchment (z ic ), and the proportion of internal catchment at the soil surface (p ic_0 ) are critical parameters in all the outputs analyzed. For imidacloprid and compound I, the most relevant parameters for drainage output are the Freundlich sorption exponent (F exp ) and z ic ; for groundwater output, the relevant parameters are F exp , the bottom depth of static macropores (z st ), and p ic_0 ; and for soil at 100-cm depth, the relevant parameters are F exp , z ic , and p ic_0 . The Morris screening and Sobol-Jansen methods match perfectly in the first position of the ranking and generally match well in the second position.
Our research results indicate that sorption, degradation, and macropore parameters related to the variation of the total relative macroporosity over depth are critical for describing pesticide leaching in macroporous soils. Therefore, the acquisition of data on pesticide and macropore parameters should be given priority when extending national soil databases with spatial soil data. Such data could improve the assessment of pesticide leaching risks in macroporous soils at regional or national scales.

C O N F L I C T O F I N T E R E S T
The authors declare no conflict of interest.