Original ResearchFrank Spurlock*
Jirka Šimůnek
Bruce Johnson
Atac Tuli

We conducted a comprehensiveMonte Carlo sensi vity analysis of
fumigant transport and vola liza on
using HYDRUS. The sensitivities of
fumigant flux and soil-gas concentraons varied by applica on scenario
and are interpreted mechanis cally.
These data provide a reference fordesigning fumigation field studies
and interpre ng their results.
F. Spurlock, B. Johnson, and A. Tuli, California
Dep. of Pes cide Regula on, 1001 I Street,
Sacramento, CA 95812; and J. Šimůnek, Dep.
of Environmental Sciences, Univ. of California–
Riverside, Riverside, CA 92521. *Corresponding
author (fcspurlock@cdpr.ca.gov).
Vadose Zone J.
doi:10.2136/vzj2012.0130
Received 7 Sept. 2012.

© Soil Science Society of America
5585 Guilford Rd., Madison, WI 53711 USA.

All rights reserved. No part of this periodical may
be reproduced or transmi ed in any form or by anymeans, electronic or mechanical, including photocopying, recording, or any informa on storage
and retrieval system, without permission in wri ng
from the publisher.

Sensi vity Analysis of Soil
Fumigant Transport andVola liza on to the Atmosphere
Interest in the use of vadose zone transport models for fumigant risk assessment is increasing. Good modeling prac ce includes an assessment of model sensi vity and output uncertainty. This computa onal study evaluated the sensi vity of HYDRUS-1D- and HYDRUS
2D/3D-simulated fumigant cumula ve flux, maximum 6-h period mean flux density, andsoil gas concentra ons to 15 model input variables using Monte Carlo La n hypercube
analysis. The input variables included fumigant physicochemical proper es, agriculturalfilm (tarp) proper es, and soil proper es. Three diﬀerent applica on scenarios were invesgated: tarped broadcast, tarped bedded shank injec on, and a tarped drip line-source
applica on. Model sensi vity to ini al water content (θ i), saturated water content (θ s), andtarp permeability varied among scenarios depending on the rela ve importance of soil
gas diﬀusive resistance and tarp mass transfer resistance to fumigant vola liza on. Model
outputs were sensi ve to fumigant air–water par on and gas-phase diﬀusion coeﬃcients,
two parameters that probably have a small contribu on to overall modeling uncertainty
because accurate es ma on methods for these parameters are available. Sensi vi es tothe fumigant degrada on rate were high in all scenarios, and sensi vity to tarp permeability was high only when substan al vola liza on occurred from the tarped por on ofsoil surfaces. Exis ng literature data for both degrada on and tarp permeability are highlyvariable; parameterizing these processes using literature es mates may contribute substan ally to model uncertainty. In several cases, the highest output sensi vi es were toθ s. For model comparisons to site-specific field data, soil texture-based es mates of θ s arepoten ally large contributors to model uncertainty; direct measurement is recommended.

Soil fumigants accounted for >20% of the 65 million kg of pesticide active
ingredients applied in 2009 California production agriculture. Application rates are on
the order of 100 kg ha−1 and methods vary, including broadcast subsurface shank injection
at 25- to 60-cm depths, bedded shank injection, and chemigation applications throughsurface and subsurface drip irrigation systems. Agricultural plastic fi lms (“tarps”) areusually used to cover the soil surface, either partially or completely, to reduce fumigant
volatilization into the atmosphere.
Volatile organic compound emissions, including those from fumigants, are regulated inareas of California that do not meet ozone air quality standards because they are ozone
precursors (Marty et al., 2010). Post-application emissions vary by fumigant and application method, and are quantified by the cumulative flux or emission ratio (cumulative flux/
fumigant applied). Peak emissions are used to estimate off-site air concentrations, which inturn serve as the basis for buffer zones to limit bystander exposure (Sullivan et al., 2004).
Agronomists concerned with efficacy typically evaluate soil gas concentrations and their
persistence in the soil profi le (Ha et al., 2009).
Field studies to estimate fumigant flux are expensive and generally yield uncertain flux estimates (Majewski, 1997; Ross et al., 1996; Wang et al., 2006; Sullivan et al., 2004). Interestin fumigant transport modeling is increasing, and several models have been used to estimate laboratory- and field-scale fumigant transport and volatilization from soils, including
CHAIN-2D (Šimůnek and van Genuchten, 1994), HYDRUS-1D (Šimůnek et al., 2008a),
HYDRUS 2D/3D (Šimůnek et al., 2011), and HWC-MODEL (Ha et al., 2009). The governing equations and associated boundary conditions for simultaneous solution of transientwater flow and convective–dispersive transport of heat and solute are quite similar among
most vadose zone transport models. Critical aspects of simulating fumigant transport in
the vadose zone, however, also include the strong temperature dependence of tarp permeabilities (Papiernik et al., 2011), post-application tarp perforation and subsequent removal,

www.VadoseZoneJournal.org

and the requirement for two volatilization boundary conditions
in situations where only a portion of the field is covered by a tarp.
Both HYDRUS-1D (Šimůnek et al., 2008a) and HYDRUS 2D/3D(Šimůnek et al., 2011) simulate tarp dynamics, while HYDRUS
2D/3D also simulates dual surface volatilization boundary conditions and complex field geometries in two or three dimensions. Several recent studies have used the HYDRUS models or their legacy
DOS-based precursor CHAIN-2D to describe fumigant transport(Luo et al., 2012; Yates et al., 2012; Luo et al., 2011; Yates, 2009;
Cryer and van Wesenbeeck, 2009, 2010).
Sensitivity analysis is a basic component of model evaluation(Crout et al., 2008; Warren-Hicks et al., 2002; van den Berg etal., 2008). Ha et al. (2009) evaluated the sensitivity of 1- and 2-dmethyl isothiocyanate soil gas concentrations to temperaturedependent partitioning and degradation parameters in a study oftarped raised beds. The degradation rate was more important thanHenry’s Law constant or the gas-phase diff usion coefficient in thatstudy. Cryer and van Wesenbeeck (2010) evaluated the sensitivityof cumulative flux to various parameter groups using the legacy
model CHAIN-2D for a single application scenario. They concluded that soil water content and tarp properties were among themost important parameters controlling cumulative flux based on
rank correlations between inputs and outputs. Neither study provided detailed sensitivities of multiple output variables in multiple
application scenarios.
The objectives of this study were to: (i) determine quantitative
sensitivities of HYDRUS-1D- and HYDRUS 2D/3D-simulated
fumigant cumulative flux, period mean flux densities, and soil gas
concentrations to 15 key model inputs in three application scenarios, (ii) interpret those sensitivities from a mechanistic standpoint,and (iii) qualitatively evaluate potential sources of uncertainty inthe model output based on the sensitivity analysis. These results
will inform sampling strategies for designed studies to compare
field-based and model-simulated flux estimates and facilitate interpretation of simulation results.

6 Methods
Overview

Two analyses were conducted in this computational study. A preliminary analysis evaluated the sensitivity of the fumigant effective
gas-phase diffusion coefficient in soil (Deff, Eq. [6] below) and theeffective soil–air surface boundary layer mass transfer coefficient(KMTC, Eq. [8] below) to seven variables (Table 1) using a Monte
Carlo Latin hypercube sensitivity analysis method. This preliminary
analysis aided interpretation of the relative importance of diffusive
resistance vs. surface mass transfer resistance to post-application fumigant volatilization in the second full sensitivity analysis.
The second analysis used the Latin hypercube method to evaluate the relative sensitivity of 50-d cumulative flux, maximum 6-h

www.VadoseZoneJournal.org

Table 1. Input variables used in the sensitivity analyses.
Description

MedianRange(min.–max.)d, cm†

boundary layer thickness425

50–800dEa , J mol−1 K−1boundary layer
activation energy−40,000−60,000 to−20,000

ΔT, °C

daily temperatureamplitude

12.5

5–20

Dg , cm2 d−1†gas-phase diff usion
coefficient

8,250

6,500–10,000

Dg Ea , J mol−1 K−1

Dg activation energy

4,650

4,500–4,800

λw, cmlongitudinal dispersivity

10.5

1–20

D w, cm2 d−1

aqueous-phase diff usion
coefficient

0.85

0.7–1.0

k1, d−1

degradation rate constant

0.36

0.03–0.69

k1 Ea , J mol−1 K−1

k1 activation energy

52,50040,000–65,000Kh†

air–water partition
coefficient

0.15

0.05–0.25

Kh Ea , J mol−1 K−1

Kh activation energy

30,00020,000–40,000

Kd, mL g−1†

soil–water partition
coefficient

0.206

0.0375–0.375

ρ b, g cm−3†

soil bulk density

1.54

1.28–1.80

θ i, cm3 cm−3†

initial soil water content

0.13

0.06–0.20

θ s, cm3 cm−3†

saturated soil watercontent

0.42

0.30–0.54

Variable†† Included in both preliminary sensitivity analyses of effective diff usion coefficient Deff and surface mass transfer coefficient KMTC and main sensitivity analysis. All others used in main sensitivity analysis only.

average flux density (μg cm2 d−1), and maximum shallow soil gas
concentration (μg cm−3) to 15 input variables. The input variables
included fumigant physicochemical properties, properties of thetarp, soil properties, and amplitude of soil surface diurnal temperature fluctuations (Table 1). Each application was simulated
to occur at 0600 h. Maximum 6-h mean flux densities were determined for post-application periods of 0600 to 1200 h, 1200 to1800 h, 1800 to 2400 h, etc., such as would be reported for a typicalfield flux experiment. Three fumigant application scenarios wereconsidered: a fully tarped broadcast application simulated with a
one-dimensional modeling domain (broadcast scenario, Fig. 1a),a tarped bedded shank application with untarped furrows simulated with a two-dimensional modeling domain (bed scenario, Fig.
1b), and an under-tarp bedded drip application simulated using aline source representation in a two-dimensional modeling domain(drip scenario) (Fig. 2).

Simula on of Water, Solute, andHeat Transport
Both HYDRUS programs simultaneously solve the Richards equation for variably saturated water flow and the advection–dispersion
equations for heat and solute transport. Solutions to those flow
and transport equations are solved numerically subject to specified

p. 2 of 12θ( h ) = θ r +θ s −θ r(1+ αh

n 1−1/n

)

;h <0
[1]

;h ≥0

= θs

where θ r is the residual water content (cm3 water cm−3 bulk soil),θ s is the saturated water content, and α (cm−1) and n are empiricalconstants (van Genuchten, 1980).
The Water Linear Reduction gas-phase tortuosity model (Moldrup
et al., 2000) was used to describe gas-phase tortuosity, τg , in all
simulations:Fig. 1. (a) One-dimensional tarped broadcast and (b) two-dimensional tarped bedded shank application modeling domains. The onedimensional figure shows nodal density on the left, initial fumigant
distribution on the right. Width (a) and diameter (b) of initial fumigant distribution was ?15 cm.

a1.5
τg = vθs[2]where θ s in this case represents total porosity and av is the air-fi lled
porosity (cm3 air cm−3 bulk soil). Tortuosity reduces the effective gas-phase diff usion of a fumigant in soil relative to that in air.
Fumigant gas-phase diffusive flux density in the soil is proportional
to the spatial concentration gradient:J =−a v τ g D g

∂C g

[3]

∂x

where J is flux density (μg cm−2 d−1) in the x direction, Dg is the
fumigant diff usion coefficient in air (cm2 d−1), and Cg is the fumigant gas-phase concentration (μg cm−3 air), which is related to thetotal soil fumigant concentration C T (aqueous + sorbed + gas, μg
cm−3 soil) through an equilibrium partition coefficient Rg:

Fig. 2. Two-dimensional tarped bed drip domain (line source
representation). Initial soil profile temperature conditions at 0600 h
on day of application shown. Defined mesh line used to evaluate water
fluxes out of the bed.

boundary conditions using Galerkin-type linear finite elements
in one (HYDRUS-1D), two, or three dimensions (HYDRUS
2D/3D) depending on the problem (Šimůnek et al., 2008b). The
heat transport equation considers conduction and convection withflowing water, while the solute transport equations consider advective and diff usive–dispersive transport in the liquid phase and diffusion in the gas phase.
For all simulations, the relationship between soil water content (θ)
and pressure head (h) was

www.VadoseZoneJournal.org

Cg =

CT
Rg

[4]

Rg =

ρbK d
θ
+
+av
Kh
Kh

[5]

where Kd (mL g−1) is the linear equilibrium soil water distribution
coefficient, ρ b is soil bulk density (g cm−3), and Kh is the dimensionless air–water partition coefficient. When the spatial gradientin C T is used in Eq. [3], the corresponding effective fumigant gasphase diff usion coefficient Deff in bulk soil is obtained by combining Eq. [2–5]:

K h D g ( θ s −θ )2.5

D eff =

θ s ⎡⎣ ρ b K d +θ + ( θ s −θ )K h ⎤⎦[6]

where θ is soil water content and av is assumed equal to θ s − θ. The
effective diff usion coefficient Deff is applicable to total fumigant
concentration and accounts for tortuosity of the air-fi lled soil porespace, the reduction in cross-sectional area due to the presence ofp. 3 of 12the soil solid phase, and fumigant partitioning among the solid,aqueous, and gas phases.

=

Boundary Condi ons
Atmospheric water flow boundary conditions (BCs) were applied at
the top surface of all modeling domains. The potential flux in that
BC is controlled by external specified time-variable precipitationand potential evaporation fluxes (Šimůnek et al., 2011). No precipitation was simulated in any of the scenarios. The potential evaporative flux was set to zero for the tarped portion of each modelingdomain, while daily potential evaporation of 6.8 mm was used forthe untarped portion of the soil surfaces. This potential evaporation
corresponds to early summer reference evapotranspiration rates in
California’s San Joaquin Valley and was simulated between 0700
and 1700 h daily, with the maximum rate between 1000 and 1600
h. For atmospheric BCs, actual evaporation is also determined by alimiting critical pressure head at the soil surface. Below that criticalpressure head, actual evaporation is less than potential evaporation.
That critical pressure head was defined as the HYDRUS default of−15,000 cm based on preliminary simulations that showed minimal
effect of that parameter. At the lower boundary, a unit hydraulicgradient (free drainage) BC was applied in all scenarios.
A stagnant layer volatile solute BC was applied at the top of all modeling domains. That BC describes fumigant surface volatilization flux
density J (μg cm−2 d−1) through a boundary of thickness d (cm) as

k MTC
Rg
K hD g[8]d ⎡⎣ ρ b K d +θ + ( θ s −θ )K h ⎤⎦Soil surface temperature was assumed to vary sinusodially with time:⎡ ⎛7 ⎞⎤T (t ) = T0 +Δ T sin ⎢ π⎜⎜⎜ 2t − ⎟⎟⎟⎥⎢⎣ ⎝12 ⎠⎥⎦[9]

where t is time (d), T0 is the average temperature, assumed here to
be 25°C, and ΔT is the amplitude of the daily fluctuation. Default
volumetric heat capacities for soil constituents and default thermal
conductivity parameters for sands (Chung and Horton, 1987) wereused. The volume fraction of soil mineral components for eachsimulation was specified as (1 − θ s).
The temperature dependence of the diffusion coefficients, air–water
distribution constant, degradation rate constant, and boundary
layer depth d were simulated using an Arrhenius-type relationship(Šimůnek et al., 2011):⎡ E ⎛ 1 1 ⎞⎤x i ,T = x i ,r exp ⎢⎢ a ⎜⎜ − ⎟⎟⎟⎥⎥⎟⎜⎣⎢ R ⎝ T r T ⎠⎦⎥

[10]

Dg

⎡ C ( 0 ) −C ( d ) ⎤
g
g
⎥⎦
d ⎢⎣
= k MTC ⎢⎡C g ( 0 )−C g (d )⎤⎥⎣⎦J=

K MTC =

[7]

where Cg(0) and Cg(d) are the solute gas-phase concentrations at
the soil surface and at the top of the equivalent boundary layerof thickness d, respectively, and kMTC (cm d−1) is a fi rst-ordermass transfer coefficient (Jury et al., 1983). The value of Cg (d)was assumed equal to zero in all cases. The equivalent boundary
layer thickness d has no physical meaning but represents a masstransfer resistance at the soil surface. Appropriate d values may
be estimated for any fumigant–tarp combination from standard
tarp kMTC measurements (Papiernik et al., 2011). For untarped
bare soil areas, d = 0.5 cm was assumed (Jury et al., 1983), while dvalues ranging up to 800 cm were used to simulate mass transfer
resistance in tarped regions based on polyethylene tarp fumigant
permeability data (Papiernik et al., 2011). At the bottom of eachmodeling domain, the solute BC was a Cauchy (third-type) BC,
where fumigant flux out of the domain is the product of the fumigant aqueous concentration and water flux.
Combining Eq. [5] and [7] gives the effective mass transfer coefficient KMTC applicable to the total fumigant concentration C T(Jury et al., 1983):

www.VadoseZoneJournal.org

where xi,T is the value of the desired coefficient at temperature T(K), xi,r is the value at the reference temperature Tr = 293 K, R is
the universal gas constant = 8.314 J mol−1 K−1, and Ea (the “activation energy,” J mol−1) is the proportionality constant between
ln(xi) and the reciprocal temperature. While the various Ea values
in Table 1 are not true activation energies in a chemical thermodynamic sense, the “activation energy” terminology is retained here
to maintain consistency with the historical Arrhenius relationship.

Fumigant Applica on Scenarios andIni al Condi ons
The geometric configuration of the one-dimensional broadcast,
two-dimensional bed, and two-dimensional drip applications,discretization, and initial fumigant concentration conditions areshown in Fig. 1 and 2. Soil profiles were uniform, and the verticalmodeling domain depths were arbitrarily chosen as 100 cm basedon preliminary simulations that demonstrated little effect of depth
on fumigant volatilization flux beyond that value. The initial profile
soil temperature in each scenario was taken as the final temperature
distribution obtained in 5.25-d preliminary simulations using an
average soil temperature of 25°C and an amplitude of 12.5°C (Eq.
[9]). In all scenarios, 100 kg ha−1 applications were simulated at
0600 h on Day 1. All simulations were conducted for an arbitrary“long” period of 50 d to ensure complete fumigant volatilization
or degradation. Observation nodes were defined in each modelingp. 4 of 12domain to obtain simulated soil gas concentrations at the end ofeach 6-h interval (Fig. 1 and 2). The observation nodes were situatedapproximately 50 cm away from the region where the fumigant wasapplied (Fig. 1 and 2). For the bed and broadcast applications, thesenodes were near the tarp surface; in the drip scenario, the nodes
were below the drip line source. The end-of-simulation fumigant
mass balance errors for 1200 simulations in each scenario, expressedas percentage of the initial application, were 0.2, 0.4, and 0.7%
(broadcast), 0.7, 1.0, and 1.5% (bed), and 0.7, 1.0, and 1.5% (drip)
for the 10th, 50th, and 90th percentile, respectively.
A majority of California fumigant applications occur in sandy
loam soils (Johnson and Spurlock, 2009), and average values forthe sandy loam texture for residual water content θ r (0.039), α(0.029 cm−1), n (1.45) (Schaap et al., 2001), and saturated hydraulic conductivity Ks (106 cm d−1) (Carsel and Parrish, 1988) wereused in all simulations. Saturated water content θ s was varied todetermine the model output sensitivity to that variable (Table 1).
In the drip scenario, the fumigant chemigation was simulated witha water application of 288 cm2 during 13.2 h, corresponding to a
drip-line water application rate of 4.4 L m−1 h−1. The application
was simulated through a 2.5-cm-diameter line source 1 cm below
the soil surface (Fig. 2). The fumigant concentration was 373 μg
L−1 and the water temperature was 12°C.

La n Hypercube Sensi vity Analysis
Latin hypercube sampling is efficient in obtaining representative
samples across the entire input parameter space (McKay et al., 1979).
The Latin hypercube sensitivity analysis method used in this study
was similar to that reported by van Griensven et al. (2006), originallyadapted from the basic method of Morris (1991). For a model output
M that is a function of N input parameters (x1, ..., xN), a sensitivity
index Si is defined that describes the effect of a perturbation Δxi ofmodel input xi, with the remaining parameters xji constant:M ( x1 ,... x i +Δ x i ,... x N )− M ( x1 ,... x i ,... x N )⎡ M ( x1 ,... x i +Δ x i ,... x N )+ M ( x1 ,... x i ,... x N )⎤ 2
⎣
⎦
xi
×
Δ xi

Si =

[11]

The individual Si are point estimates of the normalized partial
derivative of M with respect to xi, so are local measures of sensitivity. Latin hypercube sampling was used to estimate the mean
Si (denoted here as μ i) across the entire N-dimensional parameter
hyperspace. The standard deviation of Si (= σi) across that same
parameter space provides a measure of the interaction between xi
and the other inputs and of the nonlinear effect of xi on the dependent variable M (Morris, 1991; Saltelli et al., 2005).
The procedure for estimating μ i and σi was to (i) randomly select asingle value from each decile of each of the N assumed uniform xi

www.VadoseZoneJournal.org

distributions, (ii) form 10 vectors of input parameters (x1, ..., xN) by
randomly selecting one xi without replacement from each sample
drawn in the first step, (iii) run the model using each of those 10 (x1,..., xN) input vectors, and (iv) again run the model using each of the
10 input vectors but with the perturbed value xi + Δxi substituted
for xi. Thus, for each xi, 20 model runs were conducted, yielding 10
estimates of Si. The perturbations Δxi were taken as ±10% of thebase xi, with the sign of the perturbation randomly selected.
For both the preliminary and full sensitivity analyses, a single trial
consisted of estimating 10 Si for each of the input variables (Table1). Results from four trials were then used to calculate the meansensitivities μ i and standard deviation σi for each input variable–
dependent output variable combination. In the preliminary sensitivity analysis, Deff (Eq. [6], six input variables) and KMTC (Eq. [8],seven input variables) were the dependent variables (M, Eq. [11]),and that analysis was conducted using a spreadsheet.
In the second main sensitivity analysis, there were 15 input variables × 20 simulations per variable = 300 HYDRUS simulations
per trial for each modeling scenario (broadcast, bed, and drip). The
dependent variables were cumulative flux, maximum 6-h flux density, and soil gas concentration.

Input Variables Evaluated in
Sensi vity Analyses
The preliminary analysis estimated sensitivities Si of Deff (Eq. [6])and KMTC (Eq. [8]) with respect to Dg , Kh, ρ b, Kd, d, θ, and θ s(Table 1), while the second analysis evaluated the sensitivity of
HYDRUS outputs to all 15 variables listed in Table 1. Each chemical property input was assigned a uniform distribution based on theapproximate ranges of published or estimated data for 1,3-dichloropropene, chloropicrin (trichloronitromethane), methyl bromide,and methyl iodide (Table 1). Ranges for aqueous and gas-phase
fumigant diff usion coefficients (Table 1) were based on estimatesfrom the on-line SPARC calculator (http://archemcalc.com/sparc.
php, Hilal et al., 2003a, 2003b), while the gas-phase diff usion
activation energy was calculated from SPARC Dg estimates as afunction of temperature. The range of dimensionless Henry’s Law
constants were obtained from the literature (Kim et al., 2003; Glew
and Moelwyn-Hughes, 1953; Ruzo, 2006). Henry’s Law activationenergies were measured or estimated from enthalpy of vaporization data (USEPA. 2001; Kim et al., 2003; Chickos and Acree,
2003; Glew and Moelwyn-Hughes, 1953). The range of fumigant
soil degradation rates were based on the compilation of Dunganand Yates (2003), and the range of degradation activation energies
were based on the variation in aerobic degradation rates with temperature reported by Dungan and Yates (2003), Gan et al. (2000),and Guo and Gao (2009). The range for soil water partition coefficients was calculated using organic C normalized soil partition
KOC data from the California Department of Pesticide Regulation
pesticide chemistry database of registrant data submissions andthe EU FOOTPRINT database (http://sitem.herts.ac.uk/aeru/

p. 5 of 12footprint/en/), assuming a soil mass fraction organic C of 0.005.
While enhanced vapor sorption coefficients have been reported forvery dry soil conditions, we assumed this mechanism to be of minor
importance in our scenarios due to (i) fumigant label preapplication
requirements that specify soil water contents at the 20-cm depthof at least 50% of the available water content, (ii) the presence of atarp fully or partially covering the soil surface, which minimizes
evaporation, (iii) low clay contents of sandy loam soils (<20%), and(iv) the small and nonpolar nature of fumigant molecules, especiallymethyl bromide, 1,3-dichloropropene, and methyl iodide. These
factors reduce sorption enhancement effects (Petersen et al., 1996).
The uniform initial soil water content θ i spanned the approximate
range of 50 to 75% of sandy loam water contents at field capacity(h = −300 cm) as calculated from data in the UNSODA database (Leij et al., 1996). The range of saturated water contents andbulk densities were chosen from UNSODA data for sandy loam
soils. In the preliminary analysis, the uniform distribution of θ i(Table 1) was used as the sampling distribution for θ in Eq. [6]and [8]. The dispersivity maximum and minimum were based on
approximate ranges for soil columns and field soils (Radcliffe andŠimůnek, 2010, p. 249–346). The boundary layer depth and activation energy d and dEa were calculated from the polyethylene tarp
permeability data of Papiernik et al. (2011). The upper bound forthe daily temperature amplitude of 20°C was based on reported
under-tarp field temperature data (Yates et al., 2002).

6 Results

Ini al Sensi vity Analysis of theGas-Phase Diﬀusion and Boundary Layer
Mass Transfer Coeﬃcients

The sensitivity analysis of D eff (Eq. [6]) and K MTC (Eq. [8])assumed isothermal conditions, so the resulting μ i are only a general guide to the relative importance of input variables to the simulated processes. As expected from Eq. [6] and [8], the respective
sensitivities of Deff and KMTC to Kh, Dg , Kd, and ρ b were identical (Table 2). Sensitivities to ρ b and Kd were essentially identicalbecause these variables only occur together as a product in both Eq.
[6] and Eq. [8]. The principal differences between Deff and K MTC
were their sensitivities to θ and θ s, and to d, which only appears in
Eq. [8]. Both θ and θ s influence fumigant partitioning through thedenominator of both equations, but partitioning effects on Deff areovershadowed by their role in determining tortuosity. Hence, diffusive transport is particularly sensitive to both θ and θ s (Table 2).

Table 2. Mean (μ i) and standard deviation (σi) of effective gas-phasediff usion coefficient (Deff, Eq. [6]) and boundary layer mass-transfer
coefficient (KMTC , Eq.[8]) sensitivities to independent variables in the
preliminary analysis.
Deff

KMTC

Variable†

μi

σi

μi

σi

d

–

–

−0.99

0.05

ρb

−0.62

0.17

−0.62

0.17

Kd

−0.62

0.15

−0.62

0.15

θi

−1.44

0.59

−0.24

0.10

θs

2.62

0.79

−0.14

0.07

Kh

0.90

0.09

0.90

0.09

Dg

0.99

0.05

0.99

0.05

† d, boundary layer thickness; ρ b, soil bulk density; Kd, soil–water partition
coefficient; θ i, initial soil water content; θ s, saturated soil water content;
Kh, air–water partition coefficient; Dg , gas-phase diff usion coefficient.

volatilization occurred from the untarped portion of the soil surface; the median fraction of total cumulative volatilization fromthe untarped soil area was 0.90 (0.83 and 0.94 for the 10th and90th percentiles, respectively). In contrast, the fraction of total
volatilization from untarped soil in the drip scenario was only
0.11 (0.02 and 0.28 for the 10th and 90th percentiles, respectively).
Maximum 6-h mean flux densities displayed the same trend among
scenarios as the emission ratios, with medians of 56 μg cm−2 d−1
for the broadcast scenario, 132 μg cm−2 d−1 in the bed scenario,
and 286 μg cm−2 d−1 in the drip scenario (Fig. 3b). Maximum
soil gas concentrations were similar among scenarios, with medianconcentrations of 1 to 2 μg cm−3 (Fig. 3c).

General Comparison of theModeling Scenarios

General differences in temporal flux dynamics and gas concentrations at the observation points among the three application scenarios are illustrated in Fig. 4. These simulations used the medianof each input variable (Table 1). For each scenario, the temporalflux patterns across all 1200 simulations were generally similar to
those shown in Fig. 4a. In the drip scenario, the maximum 6-hmean flux density always occurred in the afternoon (1200–1800h) or evening (1800–2400 h) of the day of application, and those
fluxes were generally quite high compared with the bed and broadcast scenarios. In the bed scenario, 25% of maximum period mean
flux densities occurred during the second (1200–1800 h) period on
the day of application, but the majority (60%) were during 0600 to
1200 or 1200 to 1800 h during Day 2. The timing of the maximum
period mean flux density in the broadcast scenario was similar tothe bed scenario, but the magnitude of those peak fluxes was generally lower (e.g., Fig. 4a).

Across all 1200 simulations for each scenario, median 50-d
emission ratios for the broadcast, bed, and drip application scenarios were 0.10 (0.02 and 0.33 for the 10th and 90th percentiles, respectively), 0.21 (0.05 and 0.50) and 0.32 (0.17 and 0.61),
respectively (Fig. 3a). In the bed scenario, the majority of fumigant

Gas concentrations at the observation points were “instantaneous”
concentrations at the end of each 6-h period. While the drip scenario observation node was located 45 cm below the drip line,
these concentrations were generally comparable in magnitude to

www.VadoseZoneJournal.org

p. 6 of 12

Fig. 3. Box plots of (a) emission ratio, (b) maximum 6-h flux density,and (c) maximum fumigant soil gas concentration Cg for broadcast,bed, and drip application scenarios, showing median (center line),quartiles (box boundaries), 10th and 90th percentile (whiskers), and5th and 95th percentiles (symbols) based on N = 1200 simulations
for each scenario.

the shallow concentrations simulated in the other two scenarios(Fig. 4b). There was a strong diurnal signal in the shallow gas concentrations but much less so in the drip scenario. This was attributable to the lower temperature fluctuations at the deeper depth (Fig.
4c). The simulated concentrations were comparable to shallow gas
concentrations measured in field studies when normalized by theapplication rate (Gao et al., 2008; Gao and Trout, 2007).

www.VadoseZoneJournal.org

Fig. 4. Illustrative comparison of (a) simulated maximum 6-h mean flux
densities, (b) maximum post-application soil gas concentrations Cg , and(c) soil temperature fluctuations at the 5- and 45-cm depths. Medianvalues of input variables were used in all simulations (Table 1).

Water flow out of the bottom of the profi le was generally very low
in the broadcast and bed scenarios. In simulations using median
input variables, the cumulative water flux at 50 d was 0.035 and0.038 cm, respectively, the latter being the average flux across thebottom of the bed domain boundary. Deep drainage was higher inthe drip scenario, with a cumulative drainage flux of 0.15 cm (average across the bottom domain boundary), and this corresponded top. 7 of 12

?4% of the applied water. The maximum drainage in all scenarios
occurred at the highest θ i (0.20, Table 1). Simulations using θ i =0.20 along with the other median input variables yielded mean
cumulative water fluxes of 1.3, 1.1, and 2.0 cm for the broadcast,bed, and drip scenarios, respectively. Th is corresponds to 53% of
applied water in the drip scenario.

Sensi vity of Cumula ve Flux
While cumulative flux was relatively sensitive to the gas-phase diffusion coefficient Dg and air–water distribution coefficient Kh as
expected, the greatest sensitivity was to θ s for two of the three
application scenarios (Table 3). The μ θs for the broadcast and bed
cumulative flux (1.6 and 2.1, respectively) gives a mean change incumulative flux of 16 and 21%, respectively, for a 10% change inθ s, all other factors being equal (Table 3). In contrast, 10% perturbations to either Dg or Kh yielded nearly equivalent percentagechanges in cumulative flux for the broadcast and bed scenarios
as seen from their respective μ i values (Table 3). The high sensitivity of cumulative flux to θ s in the broadcast and bed scenarios
reflects the importance of diff usive transport in those scenarios.
The importance of diff usion is also indicated by the large negative
sensitivity of cumulative flux to θ i, consistent with preliminary
sensitivity results for Deff and KMTC (Table 2).
Sensitivities for the drip scenario were generally much lower thanfor the other two scenarios. This was partially attributable to thenear-surface application of the drip-applied fumigant directly
under the tarp and the corresponding rapid volatilization fromthe tarped area. The mean fraction of eventual total volatilizationthat occurred within the first 24 h was 0.54 for the drip application,compared with 0.24 and 0.31 for the broadcast and bed scenarios,
respectively. Thus, for the drip scenario, short mean diff usion path
lengths and rapid volatilization yielded lower flux sensitivities todiff usion-related variables: Dg , Kh, θ s, ρ b, and Kd (Table 3). Initial water content in the drip scenario had only a minor effect on
cumulative flux because the post-application water content was
dominated by the drip input, similar to the conclusion of Ha etal. (2009). Cumulative flux was sensitive to the degradation rate
constant k1 in all scenarios; however, the shorter fumigant soil residence time in the drip scenario noted above contributed to lower
μ k1 relative to bed and broadcast applications.
The sensitivity to mass transfer resistance as represented by d (Table3) varied by application method and was related to the fraction of
eventual cumulative volatilization that occurred in the tarped areaof the profi le. In the bedded shank scenario, only a small portionof the volatilization occurred through the tarp. In that scenario,
cumulative flux was essentially insensitive to d (Table 3). In contrast, cumulative flux was much more sensitive to d in the drip andbroadcast scenarios, where larger fractions of the total cumulative
volatilization occurred through tarped surfaces.

www.VadoseZoneJournal.org

Table 3. Mean (μ i) and standard deviation (σi) of cumulative flux,
maximum 6-h mean flux density, and fumigant soil gas concentration
sensitivities to 15 input variables for three application scenarios.
Input
variable†

Broadcast
μi

Bed
σi

μi

Cumulative flux
θs
1.63
1.45
2.14
Dg
1.03
0.36
0.86
Kh
0.96
0.39
0.81
dEa
0.21
0.13
0.02
Dg Ea
0.05
0.02
0.01
K h Ea
0.20
0.12
0.25
ΔT
0.01
0.08
0.00
λw
0.00
0.00
0.00
Dw
0.00
0.00
0.00
k1 Ea
−0.46
0.25
−0.41
d
−0.47
0.20
−0.04
ρb
−0.63
0.26
−0.54
Kd
−0.68
0.41
−0.58
k1
−1.06
0.49
−0.90
θi
−1.23
1.18
−1.22
Maximum 6-h mean flux density
θs
2.16
1.68
2.91
Dg
1.31
0.38
1.28
Kh
1.15
0.40
1.19
dEa
0.43
0.17
0.03
Dg Ea
0.09
0.03
0.05
K h Ea
0.28
0.17
0.49
ΔT
0.28
0.12
0.23
λw
0.00
0.00
0.00
Dw
0.00
0.00
0.00
k1 Ea
−0.26
0.22
−0.25
d
−0.54
0.19
−0.06
ρb
−0.83
0.33
−0.80
Kd
−0.92
0.46
−0.85
k1
−0.56
0.45
−0.49
θi
−1.57
1.19
−1.77
Maximum fumigant soil gas concentration
θs
1.45
1.24
1.95
Dg
0.31
0.32
0.51
Kh
1.17
0.35
1.39
dEa
−0.07
0.06
−0.11
Dg Ea
0.01
0.02
−0.02
K h Ea
0.35
0.17
0.55
ΔT
0.03
0.09
0.08
λw
0.00
0.00
0.00
Dw
0.00
0.00
0.00
k1 Ea
−0.23
0.17
−0.36
d
0.20
0.08
0.20
ρb
−0.81
0.30
−0.92
Kd
−0.84
0.41
−0.98
k1
−0.55
0.40
−0.70
−1.18
1.01
−1.60
θi

Drip
σi

μi

σi

1.30
0.35
0.36
0.01
0.02
0.12
0.02
0.00
0.00
0.25
0.01
0.24
0.38
0.51
1.05

0.44
0.48
0.47
0.16
0.03
0.18
0.12
−0.06
0.00
−0.27
−0.28
−0.17
−0.18
−0.52
−0.19

0.15
0.16
0.17
0.09
0.01
0.07
0.12
0.04
0.00
0.10
0.15
0.07
0.09
0.13
0.09

1.39
0.44
0.41
0.06
0.03
0.18
0.11
0.01
0.00
0.22
0.02
0.31
0.39
0.46
1.20

0.16
0.53
0.52
0.32
0.05
0.29
0.37
−0.10
0.00
−0.09
−0.45
−0.19
−0.20
−0.09
−0.13

0.17
0.19
0.20
0.14
0.02
0.09
0.18
0.05
0.00
0.07
0.19
0.07
0.08
0.05
0.08

1.73
0.67
0.60
0.07
0.03
0.21
0.10
0.00
0.00
0.31
0.07
0.36
0.53
0.63
1.47

−0.83
0.14
1.06
−0.02
0.00
0.25
−0.05
0.55
0.00
−0.38
0.02
−1.08
−1.10
−0.64
0.14

1.02
0.24
0.26
0.01
0.01
0.06
0.07
0.33
0.00
0.29
0.02
0.45
0.61
0.31
0.34

† θ s, saturated soil water content; Dg , gas-phase diff usion coefficient; Kh, air–water partition coefficient; dEa , boundary layer activation energy; Dg Ea ,Dg activation energy; Kh Ea , Kh activation energy; ΔT, daily temperatureamplitude; λw, longitudinal dispersivity; D w, aqueous-phase diff usion coefficient; k1 Ea , degradation rate constant activation energy; d, boundary
layer thickness; ρ b, soil bulk density; Kd, soil–water partition coefficient;
k1, degradation rate constant; θ i, initial soil water content.

p. 8 of 12

Cumulative flux was relatively insensitive to parameters that control temperature-dependent partitioning and transport (Table 3).
This supports previous suggestions that simple isothermal models
may be adequate for screening purposes if diurnal flux dynamics
are unimportant (Yates et al., 2002). In all scenarios, cumulativeflux was insensitive to Dw and λw, variables that describe transport
in the aqueous phase.

Sensi vity of Maximum Six-Hour Mean
Flux Density
While the sensitivities of 6-h mean flux densities were generally
similar to those of cumulative flux in many cases, there were somedifferences. First, sensitivities to Kd and ρ b in the broadcast andbed applications were somewhat greater. The sensitivities of thesetwo variables were essentially equal because their only computational role in HYDRUS is in phase-partitioning calculations (e.g.,
Eq. [5]). In this role, they only occur together as a product, hence afi xed percentage perturbation to either Kd or ρ b has an equivalenteffect on all model outputs.
Another difference between cumulative flux and maximum 6-hperiod mean flux densities was for the variables that control thetemperature dependence of air–water partitioning and surface
mass transfer resistance, K h Ea , dEa , and ΔT. Higher K h Ea or
ΔT yielded higher Kh, enhancing both effective diff usion (Eq. [6])and the boundary layer mass transfer coefficient (Eq. [8]). Botheffects influence period mean maximum flux density to a greater
extent than cumulative flux, for which temperature effects tended
to average out. Across scenarios, the maximum period mean flux
sensitivity μ dEa reflected the relative fraction of total volatilization associated with the tarped portion of the surface discussed
above. In contrast, the 6-h mean flux density generally displayed
low sensitivity to Dg Ea; this reflects the relatively weak temperature dependence of Dg (Table 1).

Sensi vity of Soil Gas Concentra ons
Maximum shallow gas concentrations in the broadcast and bed
scenarios displayed generally similar μ i. In these scenarios, gas concentrations were relatively insensitive to variables solely related tovolatile mass transfer across the tarped boundary layer (i.e., d anddEa) but were sensitive to those involved in partitioning (Kh, Kh Ea,
Kd, and ρ b) and diff usion (θ i, θ s, and Dg). In all cases, k1 had a negative effect on maximum shallow gas concentrations.
Soil gas concentrations were generally insensitive to ΔT, largely
due to the offsetting effect of temperature-dependent factorsthat decrease the gas-phase concentration (k1 and d) and increase
fumigant gas concentrations (Kh). Temperature was important,however, in determining the 6-h period in which maximum gas
concentrations were observed. In the bed scenario, 99% of maximum gas concentrations occurred in the 0600 to 1200 and 1200
to 1800 h time periods when surface soil temperatures were highest, while 66% of the maximum gas concentrations were similarly

www.VadoseZoneJournal.org

observed during the two high-temperature periods in the broadcast
scenario. While temperature fluctuations at the 45-cm depth were
only about 10% of those at the surface (Fig. 4c), maximum gas
concentrations at that depth in the drip scenario also occurred at
the higher temperature periods. Seventy-seven percent of maximum gas concentrations occurred during the 1800 to 2400 h time
period in that scenario. This was the time of higher soil temperatures at that depth due to the phase lag between surface and subsurface soil temperatures (Fig. 4c).
For the drip application, the dominant mechanism for fumigant
transport to the 45-cm depth under the line source was convection
in moving water as opposed to diff usion in the gas phase. This wasevident from simulation results (not shown) where Dg was set tozero. In those simulations, the median of simulated maximum gas
concentrations 45 cm below the drip emitter was only 12% lessthan in corresponding simulations using identical input variables
except with nonzero Dg. The role of convective transport in thedrip scenario was also shown by the relatively high μ λw, indicating
the contribution of convective–dispersive transport as opposed togaseous diff usion. In contrast to the broadcast and bed scenarios,
maximum gas concentrations had a negative sensitivity to θ s. This
was attributable to a greater volume of water, hence fumigant,moving deep into the profi le at low θ s (Fig. 5).

Sensi vity Nonlinearity or Interac ons
In several cases, the σi were large relative to their corresponding μ i(Table 3), reflecting nonlinear effects of variables on model output
or interactions between variables. An example of nonlinearity was
the sensitivity of the maximum 6-h mean flux density to θ i in the
bed scenario. The individual S θi were significantly correlated withθ i (Spearman r = −0.74, P < 0.001). Because S θi is essentially thederivative of the model output M with respect to θ i, it is evident
that changes in θ i lead to much greater changes in the maximum
6-h mean flux density at higher values of θ i (Fig. 6). In contrast,
SKh values for maximum 6-h flux density in the drip scenario were
only weakly dependent on Kh (Spearman r = −0.31, P = 0.052) but
showed a strong interaction with d (Fig. 7).
Two general trends in nonlinearity or interaction were evident.
First, there were significant interactions between θ i and θ s in the
bed and broadcast scenarios, where volatilization was substantiallymediated by the rate of diff usive transport. The interactive effectof these variables on flux was evident from significant Spearman
correlations between S θi and θ s, ranging from 0.39 to 0.53. Similarly, the correlations between S θs and θ i ranged from 0.41 to 0.46.
The interaction is attributable to their joint influence on Deff (Eq.
[6]). A second general result was that, while model outputs in all
scenarios displayed a negative sensitivity to the fumigant degradation rate constant k1, the effects were nonlinear. In all scenarios,
the Sk1 values of all model outputs were significantly correlated to
the k1 values themselves. Spearman correlations between Sk1 andp. 9 of 12

Fig. 6. Individual sensitivities of maximum 6-h mean flux densities (S)to initial water content θi, showing nonlinear effect of θi on maximum
flux in the bed scenario.

Fig. 5. Comparison of (a) cumulative water flux past the mesh line
below dripper (mesh line location shown in Fig. 2), and (b) soil gas
concentrations Cg across the depth range of 35 to 50 cm below thedripper for saturated water contents θs of 0.34 and 0.50.

k1 ranged from −0.36 to −0.87, indicating a greater negative effecton model outputs for shorter degradation half-lives (higher k1).

6 Discussion

The μ i are estimates of the mean normalized partial derivative ofmodel output with respect to each input. While they provide ageneral indication of variable influence on model output, knowledge of a variable’s uncertainty or variability is also required for a
relative assessment of contributions to model output uncertainty.
For example, although flux predictions were very sensitive to Dg
in some scenarios (Table 3), relatively accurate D g estimation
methods are available (e.g., Hilal et al., 2003a, 2003b; Tucker and
Nelken, 1990). Consequently, the overall contribution of Dg to theuncertainty of modeled flux estimates is probably low. Similarly,
accurate solubility and vapor pressure measurements should provide relatively accurate Kh estimates, particularly at the milligram
per liter aqueous concentration ranges expected in the field (Smithand Harvey, 2007).

www.VadoseZoneJournal.org

Fig. 7. Individual sensitivities of maximum 6-h mean flux densities toair–water partition coefficient Kh (SKh) vs. equivalent boundary layer
depth d, showing interaction of Kh with d in the drip scenario.

Most outputs were very sensitive to θ s. A common modeling practice is to estimate θ s based on soil texture class (Cryer and vanWesenbeeck, 2009; Luo et al., 2012; Yates et al., 2012). Schaap etal. (1998) reported a mean sandy loam θ s of 0.389 with a standard
deviation of 0.094 (N = 337 soils), greater than the 10% perturbation in θ s used to estimate sensitivities here. Consequently, uncertainty in θ s may contribute substantially to model output uncertainty. In addition to soil-to-soil variations within a texture class,management practices such as cultivation may result in temporalvariations in θ s within the same soil. Deviations between model
predictions obtained using θ s estimates may be potentially responsible for large differences between modeled and field-estimatedp. 10 of 12

fluxes. Alternatively, agreement between model and field-based
fluxes may be observed when the model is actually not properly
simulating field processes correctly—essentially providing the correct answer for the wrong reasons. In field studies where results areto be used as the basis for modeling, accurate site-specific determination of soil hydraulic properties is strongly preferred.
Model outputs were only moderately sensitive to the equivalentboundary layer d that describes tarp resistance to mass transfer.
Tarp mass-transfer coefficients are sometimes highly variable, however, even among similar tarp–chemical combinations. Based onan estimated diff usion coefficient in air of 8899 cm2 d−1, we calculated an iodomethane mean d of 1650 cm with a standard deviationof 1450 cm for 22 samples of unused virtually impermeable fi lm(VIF) tarp based on the diff usive resistance data of Papiernik et al.
(2011). In addition, the concordance of laboratory-measured permeabilities and actual field permeabilities are poorly understood,
particularly given the potential for stretching and tears during
application and the potentially high variation of VIF tarp permeability with humidity (Papiernik et al., 2011). Thus, for applications
where a large fraction of fumigant volatilization may occur fromtarped surfaces (e.g., broadcast or drip), accurate characterizationof tarp-specific permeabilities is important. Otherwise uncertainty
in d may substantially contribute to uncertainty in model predictions, limiting comparison with field-based flux estimates.
Model outputs were sensitive to degradation rates in all scenarios.
Even when measured under well-controlled laboratory conditions,fumigant degradation rates vary substantially with soil type, moisture content, and temperature (Dungan and Yates, 2003). Fumigant degradation is similar to tarp permeability in the sense thatthe relationship between laboratory and field conditions is poorly
understood. Consequently, degradation is potentially a major contributor to simulated estimates of flux in all scenarios.
The output variable sensitivities to the various activation energies that describe temperature dependence were generally low tomoderate. For Dg Ea and Kh Ea , estimation methods applicable
to environmental temperature ranges based on accepted physicalchemistry theory are available (Smith and Harvey, 2007; Hilal etal., 2003a, 2003b; Tucker and Nelken, 1990). Obtaining reliable
estimates for k1 Ea and dEa is more difficult. Data for the temperature dependence of lumped degradation coefficients are sparse,
although calculations for the limited fumigant aerobic degradation
data that are available suggests a relatively narrow range of 40,000
to 65,000 J mol−1 for k1 Ea may be reasonable (Dungan and Yates,
2003; Gan et al., 2000; Guo and Gao, 2009). In the case of dEa ,the sensitivity of the broadcast maximum 6-h mean flux density
was moderate. Data for dEa are also sparse. Papiernik et al. (2011)
measured the temperature dependence of the diff usion resistanceacross two VIF tarps for several fumigants. Their results showed
differences of up to a factor of two. Actual dEa measurements on
the specific tarp–fumigant combination are recommended for

www.VadoseZoneJournal.org

site-specific simulations. This is particularly true in cases where alarge fraction of volatilization occurs through the tarp (broadcast
and drip scenarios) and when period mean flux densities are ofprimary interest.
Additional factors beyond the 15 variables evaluated here also
influence flux and gas concentrations. These include chemical ororganic soil amendments (McDonald et al., 2008), water applications (Gao et al., 2008; Gao and Trout, 2007), and the depth of
fumigant application (Papiernik et al., 2004). Analysis of thesefactors was beyond the scope of this study. The sensitivity analysis
provides insight into the relative importance of different transport mechanisms under different conditions, however, and will be
useful in determining how such mitigation factors can be most efficiently implemented to reduce off-site fumigant air concentrations.

6 Conclusion
Different application scenarios yielded distinctly different sensitivity results. Interpreting the sensitivities from a mechanistic
standpoint, we concluded that the μ i are fundamentally related tothe dominant transport and volatilization mechanisms in each scenario (i.e., Eq. [6] and [8]). Volatilization in the tarped broadcast
scenario was mediated by diffusive transport and tarp mass transfer
resistances in series. In that scenario, the cumulative flux sensitivity
to boundary layer thickness was greatest, and the maximum 6-hmean fluxes were most sensitive to the boundary layer activationenergy. In contrast, the bedded tarp fluxes were less sensitive totarp parameters; volatilization occurred primarily through theuntarped furrow regions. The effect of tarping was intermediatein the drip application, where the magnitude and relative contribution of volatilization through the tarped and untarped surfaces
was dependent on both the rates of diff usive transport and soil
hydraulic characteristics that influenced post-application fumigant
and water penetration into the profi le.
Based on this sensitivity analysis, uncertainty in three areas has thegreatest potential to contribute to model uncertainty. These are (i)the fumigant degradation rate, (ii) soil water properties—both initial water content and saturated water content (i.e., porosity), and(iii) tarp permeability and permeability temperature dependence.
The latter are primarily important in scenarios such as broadcast
or near-surface drip application, where most or all volatilization
occurs from a tarped soil surface.

References
Carsel, R.F., and R.S. Parrish. 1988. Developing joint probability distribu onsof soil water reten on characteris cs. Water Resour. Res. 24:755–769.
doi:10.1029/WR024i005p00755
Chickos, J.S., and W.E. Acree, Jr. 2003. Enthalpies of vaporiza on of organic andorganometallic compounds, 1880–2002. J. Phys. Chem. Ref. Data 32:519–
878. doi:10.1063/1.1529214
Chung, S.O., and R. Horton. 1987. Soil heat and water flow with a par al surface
mulch. Water Resour. Res. 23:2175–2186. doi:10.1029/WR023i012p02175

p. 11 of 12

Crout, N., T. Kokkonen, A.J. Jakeman, J.P. Norton, L.T.H. Newham, R. Anderson,
et al. 2008. Good modelling prac ce. In: A.J. Jakeman et al., editors, Environmental modelling, so ware and decision support: State of the art and new
perspec ves. Elsevier, Amsterdam.
Cryer, S.A., and I.J. van Wesenbeeck. 2009. Es ma ng field vola lity of soil fumigants using CHAIN_2D: Mi ga on methods and comparison against chloropicrin and 1,3-dichloropropene field observa ons. Environ. Monit. Assess.
15:309–318.
Cryer, S.A., and I.J. van Wesenbeeck. 2010. Coupling field observa ons, soilmodeling, and air dispersion algorithms to es mate 1,3-dichloropropene and chloropicrin flux and exposure. J. Environ. Qual. 40:1450–1461.
doi:10.2134/jeq2010.0130
Dungan, R.S., and R.S. Yates. 2003. Degrada on of fumigant pes cides: 1,3-dichloropropene, methyl isothiocyanate, chloropicrin, and methyl bromide.
Vadose Zone J. 2:279–286. doi:10.2136/vzj2003.2790
Gan, J., S.R. Yates, F.F. Ernst, and W.A. Jury. 2000. Degrada on and vola liza onof the fumigant chloropicrin a er soil treatment. J. Environ. Qual. 29:1391–
1397. doi:10.2134/jeq2000.00472425002900050004x
Gao, S., and T.J. Trout. 2007. Surface seals reduce 1,3-dichloropropene and chloropicrin emissions in field tests. J. Environ. Qual. 36:110–119. doi:10.2134/
jeq2006.0107
Gao, S., T.J. Trout, and S. Schneider. 2008. Evalua on of fumiga on and surface
seal methods on fumigant emissions in an orchard replant field. J. Environ.
Qual. 37:369–377. doi:10.2134/jeq2007.0088
Glew, D.N., and E.A. Moelwyn-Hughes. 1953. Chemical sta cs of the methyl halides in water. Discuss. Faraday Soc. 15:150–161. doi:10.1039/df9531500150
Guo, M., and S. Gao. 2009. Degrada on of methyl iodide in soil: Eﬀects of environmental factors. J. Environ. Qual. 38:513–519. doi:10.2134/jeq2008.0124
Ha, W., C.D. Stanley, R.S. Mansell, and H.A. Ajwa. 2009. 2-D Simula on of nonisothermal fate and transport of a drip-applied fumigant in plas c-mulchedsoil beds: II. Sensi vity analysis and model applica on. Transp. Porous Media 76:431–448. doi:10.1007/s11242-008-9256-2
Hilal, S.H., S.W. Karickhoﬀ, and L.A. Carreira. 2003a. Predic on of chemical reac vity parameters and physical proper es of organic compounds from molecular structure using SPARC. EPA/600/R-03/030. USEPA Natl. Expos. Res.
Lab., Research Triangle Park, NC.
Hilal, S.H., S.W. Karickhoﬀ, and L.A. Carreira. 2003b. Verifica on and valida onof the SPARC model. EPA/600/R-03/033. USEPA Natl. Expos. Res. Lab., Research Triangle Park, NC.
Johnson, B., and F. Spurlock. 2009. Dominant soil types associated with fumigant applica ons in ozone nona ainment areas. Technical memorandum.
Dep. Pest. Reg., Sacramento, CA. h p://www.cdpr.ca.gov/docs/emon/pubs/
ehapreps/eh4342.pdf.
Jury, W.A., W.F. Spencer, and W.J. Farmer. 1983. Behavior assessment modelfor trace organics in soil: I. Model descrip on. J. Environ. Qual. 12:558–564.
doi:10.2134/jeq1983.00472425001200040025x
Kim, J., S.K. Papiernik, W.J. Farmer, J. Gan, and S.R. Yates. 2003. Eﬀect of formula on on the behavior of 1,3-dichloropropene in soil. J. Environ. Qual.
32:2223–2229. doi:10.2134/jeq2003.2223
Leij, F.J., W.J. Alves, M.Th. van Genuchten, and J.R. Williams. 1996. The UNSODA unsaturated soil hydraulic database: User’s manual version 1.0.
EPA/600/R-96/095. USEPA Natl. Risk Manage. Res. Lab., Cincinna , OH.
Luo, L., D.A. Ashworth, J. Šimůnek, R. Xuan, and S.R. Yates. 2012. Assessment of methods for methyl iodide emission reduc on and pest control
using a simula on model. Atmos. Environ. 66:33–40. doi:10.1016/j.atmosenv.2012.03.080
Luo, L., S.R. Yates, and D.J. Ashworth. 2011. Predic ng methyl iodide emission,soil concentra on, and pest control in a two-dimensional chamber system. J.
Environ. Qual. 40:109–117. doi:10.2134/jeq2010.0280
Majewski, M.S. 1997. Error evalua on of methyl bromide aerodynamic flux
measurements. ACS Symp. Ser. 652:135–153.
Marty, M., F. Spurlock, and T. Barry. 2010. Vola le organic compounds frompes cide applica on and contribu on to tropospheric ozone. In: R. Krieger,
editor, Hayes’ handbook of pes cide toxicology. 3rd ed. Academic Press,
London. p. 571–586.
McDonald, J.A., S. Gao, R. Qin, T.J. Trout, and B.D. Hanson. 2008. The eﬀectof chemical and manure amendment with water applica on and tarp on
1,3-dichloropropene emissions from soil. Environ. Sci. Technol. 42:398–402.
doi:10.1021/es071133o
McKay, M.D., R.J. Bechman, and W.J. Conover. 1979. A comparison of threemethods for selec ng values of input variables in the analysis of output
from a computer code. Technometrics 21:239–245.
Moldrup, P., T. Olesen, J. Gamst, P. Schjonning, T. Yamaguchi, and D.E. Rolston.
2000. Predic ng the gas diﬀusion coeﬃcient in repacked soil: Water-induced
linear reduc on model. Soil Sci. Soc. Am. J. 64:1588–1594. doi:10.2136/
sssaj2000.6451588x
Morris, M.D. 1991. Factorial sampling plans for preliminary computa onal experiments. Technometrics 33:161–174. doi:10.1080/00401706.1991.10484804
Papiernik, S.K., R.S. Dungan, W. Zheng, M. Guo, S.M. Lesch, and S.R. Yates. 2004.
Eﬀect of applica on variables on emissions and distribu on of fumigants
applied via subsurface drip irriga on. Environ. Sci. Technol. 38:5489–5496.
doi:10.1021/es049064q

www.VadoseZoneJournal.org

Papiernik, S.K., S.R. Yates, and D.O. Chellimi. 2011. A standardized approachfor es ma ng the permeability of plas c films to soil fumigants undervarious field and environmental condi ons. J. Environ. Qual. 40:1375–1382.
doi:10.2134/jeq2010.0118
Petersen, L.W., P. Moldrup, O.H. Jacobsen, and D.E. Rolston. 1996. Rela ons
between specific surface area and soil physical and chemical proper es. Soil
Sci. 161:9–21. doi:10.1097/00010694-199601000-00003
Radcliﬀe, D.E., and J. Šimůnek. 2010. Soil physics with HYDRUS: Modeling andapplica ons. CRC Press, Boca Raton, FL.
Ross, L.J., B. Johnson, K.D. Kim, and J. Hsu. 1996. Predic on of methyl bromide
flux from area sources using the ISCST model. J. Environ. Qual. 25:885–891.
doi:10.2134/jeq1996.00472425002500040033x
Ruzo, L.O. 2006. Physical, chemical and environmental proper es of selectedchemical alterna ves for the pre-plant use of methyl bromide as soil fumigant. Pest Manage. Sci. 62:99–113. doi:10.1002/ps.1135
Saltelli, A., M. Ra o, S. Tarantola, and F. Campolongo. 2005. Sensi vity analysis
for chemical models. Chem. Rev. 105:2811–2828. doi:10.1021/cr040659d
Schaap, M.G., F.J. Leij, and M.Th. van Genuchten. 2001. ROSETTA: A computer program for es ma ng soil hydraulic parameters with hierarchical pedotransfer func ons. J. Hydrol. 251:163–176. doi:10.1016/S00221694(01)00466-8
Schaap, M.G., F.J. Leij, and M.Th. van Genuchten. 1998. Neural network analysis for hierarchical predic on of soil hydraulic proper es. Soil Sci. Soc. Am. J.
62:847–855. doi:10.2136/sssaj1998.03615995006200040001x
Šimůnek, J., M. Šejna, H. Saito, M. Sakai, and M.Th. van Genuchten. 2008a. The
HYDRUS-1D so ware package for simula ng the one-dimensional movementof water, heat, and mul ple solutes in variably-saturated media. Version 4.0.
HYDRUS So w. Ser. 3. Dep. of Environ. Sci., Univ. of California, Riverside.
Šimůnek, J., and M.Th. van Genuchten. 1994. The CHAIN-2D code for simulating two-dimensional movement of water flow, heat, and mul ple solutes
in variably-saturated porous media. Version 1.1. Res. Rep. 136. U.S. Salinity
Lab., Riverside, CA.
Šimůnek, J., M.Th. van Genuchten, and M. Šejna. 2008b. Development andapplica ons of the HYDRUS and STANMOD so ware packages, and related
codes. Vadose Zone J. 7:587–600. doi:10.2136/VZJ2007.0077
Šimůnek, J., M.Th. van Genuchten, and M. Šejna. 2011. The HYDRUS so ware
package for simula ng two- and three dimensional movement of water,
heat, and mul ple solutes in variably-saturated media. Version 2.0. PC-Progress, Prague, Czech Republic.
Smith, F.L., and A.H. Harvey. 2007. Avoid common pi alls when using Henry’sLaw. Chem. Eng. Prog. 103:33–39.
Sullivan, D.A., M.T. Holdsworth, and D.J. Hlinka. 2004. Monte Carlo-based dispersion modeling of oﬀ-gas releases from the fumigant meta-sodium for determining distances to exposure endpoints. Atmos. Environ. 38:2471–2481.
doi:10.1016/j.atmosenv.2004.02.001
Tucker, W.A., and L.H. Nelken. 1990. Diﬀusion coeﬃcients in air and water. In:
W.J. Lyman et al., editors, Handbook of chemical property es ma on methods. Am. Chem. Soc., Washington, DC. p. 17-1–17-24.
USEPA. 2001. Correc ng the Henry’s Law constant for soil temperature. Fact
Sheet. USEPA Ecosyst. Res., Athens, GA. h p://www.epa.gov/athens/learn2model/part-two/onsite/doc/factsheet.pdf.
van den Berg, F., D.J. Brus, S.L.G.E. Burgers, G.B.M. Heuvelink, J.G. Kroes, J.
Stolte, et al. 2008. Uncertainty and sensi vity analysis of GeoPEARL. Alterra
Rep. 1330. PBL Rep. 500123001. Alterra, Wageningen, the Netherlands.
h p://www.pearl.pes cidemodels.nl/pdf/report_1330.pdf
van Genuchten, M.Th. 1980. A closed-form equa on for predic ng the hydraulic conduc vity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892–898.
doi:10.2136/sssaj1980.03615995004400050002x
van Griensven, A., T. Meixner, S. Grunwald, T. Bishop, M. Diluzio, and R.
Srinivasan. 2006. A global sensi vity analysis tool for the parameters of
mul -variable catchment models. J. Hydrol. 324:10–23. doi:10.1016/j.jhydrol.2005.09.008
Wang, L., D.B. Parker, C.B. Parnell, R.E. Lacey, and B.W. Shaw. 2006. Comparison of CALPUFF and ISCST3 models for predic ng downwind odor andsource emission rates. Atmos. Environ. 40:4663–4669. doi:10.1016/j.atmosenv.2006.04.043
Warren-Hicks, W., J.P. Carbone, and P.L. Havens. 2002. Using Monte Carloanalysis to judge model predic on accuracy: Valida on of the Pes cide
Root Zone Model 3.12. Environ. Toxicol. Chem. 21:1570–1577. doi:10.1002/
etc.5620210807
Yates, S.R. 2009. Analy cal solu on describing pes cide vola liza on fromsoil aﬀected by a change in surface condi on. J. Environ. Qual. 38:259–267.
doi:10.2134/jeq2008.0059
Yates, S.R., D. Ashworth, and L. Luo. 2012. Reducing emissions of vola le organic compounds (VOCs) from agricultural soil fumiga on: Comparing emission es mates using simplified methodology. U.S Salinity Lab., Riverside, CA.
www.arb.ca.gov/research/rsc/10-28-11/item4dfr07-332.pdf
Yates, S.R., D. Wang, S.K. Papiernik, and J. Gan. 2002. Predic ng pes cide volaliza on from soils. Environmetrics 13:569–578. doi:10.1002/env.542

p. 12 of 12