Root water uptake of biofuel crops revealed by coupled electrical resistivity and soil water content measurements

Biofuel crops, including annuals such as maize (Zea mays L.), soybean [Glycine max (L.) Merr.], and canola (Brassica napus L.), as well as high‐biomass perennial grasses such as miscanthus (Miscanthus × giganteus J.M. Greef & Deuter ex Hodkinson & Renvoiz), are candidates for sustainable alternative energy sources. However, large‐scale conversion of croplands to perennial biofuel crops could have substantial impacts on regional water, nutrient, and C cycles due to the longer growing seasons and differences in rooting systems compared with most annual crops. However, due to the limited tools available to nondestructively study the spatiotemporal patterns of root water uptake in situ at field scales, these differences in crop water use are not well known. Geophysical imaging tools such as electrical resistivity (ER) reveal changes in water content in the soil profile. In this study, we demonstrate the use of a novel coupled hydrogeophysical approach with both time domain reflectometry soil water content and ER measurements to compare root water uptake and soil properties of an annual crop rotation with the perennial grass miscanthus, across three growing seasons (2009–2011) in southwest Michigan, USA. We estimated maximum root depths to be between 1.2 and 2.2 m, with the vertical distribution of roots being notably deeper in 2009 relative to 2010 and 2011, likely due to the drought conditions during that first year. Modeled cumulative ET of both crops was underestimated (2–34%) relative to estimates obtained from soil water drawdown in prior studies but was found to be greater in the perennial grass than the annual crops, despite shallower modeled rooting depths in 2010 and 2011.

high-biomass, productive perennial grasses, woody plants, or non-grain-crop residues such as maize (Zea mays L) stover that can be harvested for bioenergy production as an alternative to petroleum fuel sources. Leading candidates of perennial grasses include miscanthus (Miscanthus ×giganteus J.M. Greef & Deuter ex Hodkinson & Renvoiz) and switchgrass (Panicum virgatum L.), which have high water use efficiency and lower N demand relative to maize (Stenjem et al., 2019). Socioeconomic models have projected that >3 million ha of cropland in the United States will need to be converted to energy crop production to meet state and federal standards by 2030 (Oliver & Khanna, 2017). The environmental impacts of such large-scale land use change potentially include altered regional water balances (Abraha et al., 2015;Hickman et al., 2010), although that depends on many factors, including the soil and climatic conditions, crop varieties and cultivars, and agronomic management. There is thus a need to anticipate how an emerging large-scale shift in land use to perennial biofuel crops may affect cropland water balances, along with resulting subsurface and surface water resources.
In the face of the projected growth in biofuel crop production and its potential consequences for water resources, new tools are needed to measure and simulate how crop water use will respond to agricultural management, soil properties, and climate variability. Modeling these effects is currently limited by gaps in field data for model validation and the knowledge of the biophysical processes that drive water and nutrient cycling in these nascent agricultural systems (Uhlenbrook, 2007). Long-term field measurements of water balance for various candidate biofuel cropping systems are limited to a handful of study sites and do not cover the range of environments where these crops are likely to be grown (Hamilton et al. 2015). Studies that have quantified the likely effects on the water balance of land use conversion to biofuels in the midwestern United States have included watershed to regionalscale modeling scenarios, as well as field-scale measurements of evapotranspiration (ET). Modeling results have suggested that replacement of annual with perennial crops, such as the candidate grasses for biofuel production, will likely increase ET and would thus reduce groundwater recharge and overland runoff (Georgescu et al. 2011;Schilling et al. 2008;Vanloocke et al. 2010). Side-by-side comparison of perennial grasses and maize in Illinois found ET to be higher in the perennial grasses (Hickman et al. 2010). However, Hamilton et al. (2015) and Abraha et al. (2020) found that this was not necessarily the case at a southern Michigan site, where ET from the perennial grasses was comparable to maize in both wet and dry years based on measurements of soil water uptake and eddy covariance, respectively. They proposed that this discrepancy could be attributed to development of water-limited conditions during most years in the well-drained sandy loam soils of the Michigan site.

Core Ideas
• Expansion of cellulosic biofuels may cause large impacts on the water balance. • Root processes that control ET are poorly characterized in new agricultural systems. • Electrical resistivity measurements are sensitive to root water uptake processes. • A new coupled hydrogeophysical inversion can estimate crop parameters. • Higher LAI and earlier emergence results in greater modeled ET in the biofuel grasses.
At the plant scale, root function research has mainly been limited to annual plants, and thus much less is known about the spatiotemporal patterns of root water uptake by perennial grass crops (Ryan et al. 2016). Understanding the likely implications for water balances of conversion from annual crops or native grasslands to biofuel cropping systems requires knowledge of how root water uptake and transpiration compare between the preexisting vegetation and the new cropping system. This matter is complicated by the variability in root distribution behavior across plant species; a study by Mann et al. (2013) found that switchgrass and miscanthus had notable differences in root growth in response to drought and that the rooting depth of miscanthus increased significantly under irrigated conditions. Since growing biofuel crops on marginal soils is a desirable strategy to avoid competition with food crops on prime farmland, studying differences in root adaptations to dry conditions is of particular importance for accurately modeling ET.
Existing methods to characterize spatiotemporal variation in root biomass and distribution are impractical beyond the scale of individual plants because measuring roots in situ is labor intensive and destructive. Coring, trenches, or excavation of the root zone can provide an indication of the biomass distribution with depth (Neukirchen et al., 1999), and root cores can measure fine root growth (Smit et al., 2000). However, traditional extraction and sieving of roots may not capture all fine roots, nor can it discriminate between active and inactive roots. Root windows and photography have been used to monitor root growth of several prospective biofuel crops (Mann et al., 2013), but this approach provides only a one-or two-dimensional (1D or 2D) sampling of a restricted area.
Electrical resistivity (ER) measurements provide a novel and minimally invasive hydrogeophysical method to study root-water interactions (Cimpoiaşu et al., 2020). Resistivity of the soil-plant-water continuum varies temporally with fluctuations in soil water content, as well as temperature and fluid conductance. Recent research has taken advantage of this sensitivity to identify and quantify spatial zones of root water uptake from the decrease in soil water content and concomitant increase in ER over the growing season (Bass et al., 2017;Garré et al., 2013;Jayawickreme et al., 2008Jayawickreme et al., , 2010Robinson et al., 2012;Vanella et al., 2018). These methods all rely on the use of geophysical models to "invert" the ER data (i.e., to transform the measured apparent resistivity values into laterally and depth-variable values of resistivity, which can then be related to soil water content). More recently, coupled hydrogeophysical inversion methods using both forward hydrologic and geophysical models have been applied to study shallow subsurface water dynamics including hydrological model parameterization (Kuhl et al., 2018;Moreno et al., 2015;Tran et al., 2016). Using ER to build a process-based model of the soil water dynamics enables the user to test hypotheses regarding the drivers of the observed changes, as well as the ability to test future and control scenarios. This approach also avoids nonunique and unconstrained solutions inherent to smoothing processes in traditional ER data inversions (Hinnell et al., 2010).
The coupled hydrogeophysical inversion approach also has advantages over calibrating a hydrologic model only with measurements made by soil water content probes because it provides a multidimensional indication of changes in soil water content. Although soil water content measurements accurately represent in situ conditions, their coverage is often limited to profiles at single points in a field. Previous research has demonstrated the use of point measurements to estimate root water uptake (Hamilton et al., 2015) and calibrate hydrological models (Schelle et al., 2012), but such approaches are typically not able to characterize within-field heterogeneities.
Here, we use for the first time a coupled hydrogeophysical inversion approach to estimate ET and root depth from ER and soil water content measurements across several nonirrigated biofuel crops in southwestern Michigan, USA, over a 3-yr period from 2009 to 2011, working in the same experiment where Hamilton et al. (2015) estimated ET from soil water profiles. We develop hydrogeophysical models for two contrasting cropping systems, one of an annual crop rotation and another of the perennial miscanthus grass, to estimate parameters controlling root water uptake and soil petrophysics with a global optimization algorithm. We examine soil water heterogeneity across the field site and consider its implications for cumulative ET for each crop throughout the model period. Although our results are specific to the observed cropping systems, this approach should be broadly applicable for any plant-soil-atmosphere system with sufficient data to build a representative hydrological model and soil electrical properties appropriate for ER methods.

Study site
Our study site is the Great Lakes Bioenergy Research Center, an experimental facility on a glacial outwash plain in southwestern Michigan, USA ( Figure 1a). The region's climate is temperate and humid, with mean annual precipitation of 900 mm and mean annual temperature of 9˚C. The soils in the region are classified as Typic Hapludalfs, with a bimodal distribution of silt and sand in the upper 0.5 m due to loess deposition that holds more water than the underlying coarse sand parent material (Luehmann et al., 2016). The thick glaciofluvial deposits onsite are rich in carbonates (calcite, dolomite) below the chemical weathering front at ∼2m depth. At this interface, CO 2 driven carbonate dissolution releases Ca 2+ , Mg 2+ , and HCO 3 ions resulting in porewater electrical conductivities ranging between 400 and 700 μS cm −1 and an average groundwater electrical conductivity of 550 μS cm −1 . The water table is at a depth of ∼17 m.
Established in 2008, the Biofuel Cropping System Experiment at the research center maintains five replicate blocks of 10 different biofuel crop treatments on individual 28-m × 40m plots (Figure 1b). For this study, we investigated two of the treatments in Block 1, comparing root water dynamics of an annual crop rotation ( The annual crops were grown using conventional agronomic management for the region (Sanford et al., 2016). All crops were grown without tillage and aboveground biomass was harvested annually. Miscanthus rhizomes were planted with 0.75m spacing within and between rows for a total of 52 plants per row in 36 total rows within the plot. Canola was planted at a rate of 9.75 kg seeds ha −1 with 0.20-m row spacing, maize at a rate of 70,000 seeds ha −1 with 0.75-m row spacing, and soybean at a density of 462,500 seeds ha −1 with 0.40-m row spacing.

2.2
Climate, soil, and vegetation data Soil temperature and water content profiles were collected at a single location roughly 13 m from the ER array and roughly 5 m from the edge of each plot. Temperature was measured hourly using data-logging temperature sensors (Thermochron iButton DS1922L) installed vertically at three depths (0.24, 0.64, and 1.25 m) along a 5-cmdiam. buried polyvinyl chloride (PVC) pipe. Soil water content was recorded hourly using time domain reflectometry (TDR) probes inserted horizontally at six depths from 0.2 to 1.25 m. For a full site description and details of the TDR installation and calibration methods, see Hamilton et al. (2015). Raw hourly TDR data were filtered using a 24-h moving median and moving standard deviation to reduce data noise resulting in 212,718 datapoints (78% recovery).
Soil samples were taken from each plot to characterize texture and density of the soil profiles. Samples were collected with a bucket auger and aggregated into seven sections, from 0 to 0.1 m and then at six 0.2-m intervals to a depth of 1.3 m. Percentage sand, silt, and clay were measured within each sample using a Malvern Mastersizer laser. Bulk density (g cm −3 ) was estimated from the dry weight of each section. Estimates of field capacity and wilting point for each soil layer were visually inferred from the 2009 TDR data, which experienced sig-nificant drying throughout July and early August. Soil physical properties are shown in Table 1.
Leaf area index (LAI) values were measured for each crop in all four cardinal directions at weekly intervals from the period when aboveground growth first appeared until plant senescence in the fall during the 2010 season using a handheld Li-Cor instrument. To prepare the hourly model inputs (described below), the noisy LAI data were smoothed by taking the mean of the four directional measurements, then calculating a 5-wk moving median. A linear interpolation was then used to downscale from weekly to hourly LAI for each crop in 2010. Much less frequent LAI measurements were available for 2009 and 2011, and therefore the same interpolated LAI curve from 2010 was used to approximate 2009 and 2011 hourly LAI values, adjusting the curve forward or backward in time to account for shifts in planting day for the annual crops and emergence (determined from plant cameras recording daily) in the perennial crop. The 2010 LAI data for the canola and soybean rotations were measured at two additional treatment plots located directly adjacent to Plot G3 ( Figure 1) and were used to approximate the 2009 and 2011 hourly LAI used in the model.
Depth profiles of root biomass data were available only for the perennial miscanthus; unfortunately, only total root biomass was measured for the annual crop rotation. Root biomass for miscanthus was sampled across four depth zones (0-0.10, 0.10-0.25, 0.25-0.50 ,and 0.50-1.0 m), using three 0.06-m-diam. cores taken along a gradient from the plant: at the center of the plant bunch, at the edge of the plant, and in the space between the plants along the rows. Each core was taken after senescence at the end of each growing season in each of the five replicate miscanthus plots. A proxy for root mass density (g m −3 ) distribution with depth was calculated T A B L E 1 Summary of soil properties, including soil texture, percentage sand, silt, and clay, bulk density, field capacity, wilting point, and the corresponding soil box resistivity at 25˚C (all other soil box resistivity data are presented in Supplemental Figure S1). Note the deeper transition from loamy sand to sand at the annual rotation plot and the increase in resistivity from sandy loam to loamy sand, as well as the distinctly lower resistivity in the carbonate-rich sand vs. by dividing the total measured root mass (g m −2 ) in each core by the core length (m).

Electrical resistivity and soil petrophysical data
Prior to the 2009 growing season, graphite rod electrodes (chosen over other materials like stainless steel to limit corrosion) were permanently installed within each experimental plot. Within each plot, 40 rods, each 0.08 m long and 12 mm in diameter, were placed in a linear array with 0.3m spacing in trenches (see Figure 1b, c), perpendicular to the crop rows, with the top of the electrodes at ∼0.3-m depth (a similar approach is described in Blanchy et al., 2020). Electrodes were pushed down into undisturbed soil at the bottom of the trench and were intentionally buried below the plow zone to avoid disturbance by farming activities such as planting and harvest. The rods were wired to a takeout box outside each experimental plot. The relative locations and elevations of each electrode were surveyed using a total station prior to trench backfilling. The disturbance of the soil due to the trenching and backfilling process was assumed to have a negligible effect relative to the history of agricultural activities at the study site.
The ER data were collected from the permanent arrays with an AGI SuperSting R8 approximately monthly during the growing seasons from 4 May to 22 Sept. 2009, 3 May to 20 Sept. 2010 and 2 May to 21 Sept. 2011. Each survey included 439 measurements of resistance, R (Ω) collected in a dipole-dipole configuration, where two current electrodes, C 1 and C 2 were placed left of two potential electrodes, P 1 and P 2 at a range of electrode spacings to create a 2D pseudosection, where and ΔV (V) is the measured potential drop across electrodes, and I (A) is the applied current. Greater prospecting depths are achieved through increasingly spaced current and potential electrode pairs. The effective depth, D e (m), is approximated as 0.2 times the distance between the C 1 and P 1 electrodes, which in this study ranged from 0.6 to 9.6 m, with current pair spacing ranging from 0.3 to 1.2 m. Measurement errors caused by electrode connection failures during a survey were later removed by a two-step filtering process of values ± three standard deviations from the moving lateral median from that survey for each unique effective depth, and then repeated with ± one standard deviation. Across all surveys, the median data removal with this method was 20 measurements per survey, or around 5%. For the 1D parameter estimation portion of the study, filtered data were collapsed from 2D to 1D by extracting a median value for each unique effective depth. Soil-box resistivity tests were conducted on soil samples extracted at 0.08-m increments from two 3.6-m-long x 0.038-m-diam. cores taken with a Geoprobe MT540 at an adjacent plot with comparable soil. We followed the standard protocol for soil-box resistivity measurements established in ASTM G187. Samples were prepared for placement in the soil-box by drying the soil at 105˚C for 24 h to remove all water, sieving at 1 mm, and adding distilled water in increments equivalent to volumetric water contents of 0.05, 0.10, 0.15, 0.20, and 0.30 cm cm −3 (Table 1). Soil box resistivity measurements within each soil texture class were averaged to estimate the relationship between ER and water content at the wilting point and field capacity for each soil texture. Only one sample was available from the loamy sand soil class, therefore confidence intervals are not reported for that layer (Supplemental Figure S1).

Coupled hydrogeophysical inversion
To study the root-water dynamics of each crop type, we used a novel coupled hydrogeophysical inversion method with itera-tive optimization enhanced from that presented in (Kuhl et al., 2018). For coupled hydrogeophysical inversion, rather than inverting the static ER data to obtain static soil water content distributions, a transient hydrological model is coupled with a forward geophysical model, and model parameters are optimized to match observed geophysical and hydrological data. This inversion method offers numerous advantages, chiefly that assumptions often required for traditional inversion such as surface-placed electrodes, topographic invariance, or material lateral homogeneity are relaxed. Indeed, the complexity of the geologic material is limited only by the complexity of the forward hydrologic and geophysical models. The specific steps for the coupled inversion applied here are (a) a transient hydrological model is run across all observation times, including spin-up periods needed for model equilibration; (b) static 1D soil water content profiles are extracted from the transient hydrological model at times corresponding with the date and time of the ER surveys; (c) a layer-specific petrophysical relationship is applied to convert soil water content to first reference-temperature resistivity, then corrected for the in situ temperature; (d) the 2D electrical potential field is then calculated in a forward ER model to compare with the measured resistances from the ER survey; and lastly, (e) the hydrologic and petrophysical model parameters are then optimized to fit observed TDR and ER data using a global optimization algorithm.
In this study, we introduce a new two-step optimization routine that allows for the model to invert lateral soil-type transitions using the ER data. For this procedure, we first optimize the hydrologic, root, and petrophysical model parameters using lateral averages of both modeled and observed ER data. For this first step, we use a single 1D hydrologic model with a fixed soil texture interface depth inferred from soil texture data. Second, using those optimized parameters, we then run multiple 1D hydrologic models, each applied to an interval along the ER transect. Soil layer transitions are optimized for each 1D hydrologic model to account for the alongtransect lateral variability in soil properties.

Forward hydrological model
We built a HYDRUS-1D soil hydrological model (Šimůnek et al., 2009) for each crop to simulate the hourly 1D hydrological fluxes of each plot for a 3-yr model period from 1 Nov. 2008 to 1 Nov. 2011. HYDRUS-1D primarily solves Richards' equation for variably saturated flow with many options for simulating additional subsurface transport processes. For this study, we ran HYDRUS with root growth, root water uptake (transpiration), evaporation, heat and CO 2 respiration and transport (Suarez & Šimůnek, 1993), and snow hydrology, as well as the Unsatchem module (Suarez & Šimůnek, 1997) to model soil water ion concentrations. Major ion chemistry was incorporated into this model primarily to account for the known spatial and temporal variability in pore water ionic concentrations due to seasonal CO 2 driven carbonate dissolution around 2-m depth. The model was vertically discretized into a 0.04-m grid down to 17-m depth, with an hourly output. Upper and lower boundary conditions were specified for each of the water, heat, and CO 2 transport modules. Hourly atmospheric fluxes from the weather station, including precipitation (cm h −1 ) and potential ET(cm h −1 ) (partitioning method described below) with a surface runoff threshold of 2 mm, and free drainage, were used as the upper and lower hydrologic boundary conditions, respectively. Canopy storage (mm) was calculated as 15% of LAI (Dai et al., 2003;Dickinson et al., 1991) and was used to reduce incoming precipitation until maximum storage capacity (∼1 mm) was exceeded during individual rain events (defined as consecutive hours with measured precipitation). Hourly ground surface temperatures from the weather station were used as the upper boundary condition of the heat transport model, and a heat flux boundary condition was set for the bottom thermal boundary. Atmospheric concentrations of CO 2 (0.00033 m 3 m −3 ) were used as the upper boundary condition in the CO 2 transport model, whereas a zero gradient condition was used for the bottom boundary. Root respiration of CO 2 was modeled using the default parameters provided in HYDRUS.
Calculated potential reference grass ET, ET 0REF (cm h −1 ), was adjusted to a crop specific potential ET, ET 0C (cm h −1 ), using an estimated maximum crop coefficient K C [-] that was temporally scaled with LAI (Equation 4) and partitioned into potential evaporation, E 0 , and potential transpiration, T 0 (cm h −1 ), using the interpolated hourly LAI curves for each crop and a soil cover fraction, SCF [-] (Equations 2-5) (Šimůnek et al., 2009). During nongrowing periods, K C was set equal to 0.4 (Allen et al., 1998). In the absence of site data to the contrary, nighttime transpiration was assumed to be negligible and was set to zero when solar radiation was zero.
Model structural parameters were then specified as follows: for the soil hydraulic model, we used the van Genuchten-Mualem equation. Parameters of the hydraulic model for each layer were derived from the Rosetta database (Schaap et al., 2001) using the grain size analysis and bulk density measurements for each plot (Table 1). Initial soil layer transition depths (distinct from the HYDRUS computational layers) were set as half the distance between the depths of the observed changes in grain size unless the TDR data suggested a layer boundary should be adjusted (transition depths shown in Table 1). Heat transport parameters were calibrated with soil temperature data from previous work at the study site provided in Kuhl et al. (2018). Parameters controlling CO 2 transport and production were also taken from the database in HYDRUS but were not modified.
We used a separate HYDRUS 1-D executable including a new root growth module (Hartmann et al., 2018) that allows increased flexibility of the modeled shape and depth of the root distribution through time. We selected the Vrugt (Vrugt et al., 2001) equation to control the shape of the root density distribution, beta: Where beta [-] is the root density at depth z (m), R D (m) is the maximum rooting depth, and pz [-] and zv (m) are fitting parameters that control the shape of exponential decay in the root distribution. The root growth, considered here as the increase in root depth with time, is given by the following equation (Hao et al., 2005): where dd (m) is the potential increase in rooting depth at time t, t p (d) and t m (d) are the time to planting and time to maturity of the plant, respectively, Z (m) is the potential rooting depth from time t -1, and R D (m) is, again, the static maximum rooting depth. The t m parameter can be altered to change how quickly the plant reaches its maximum rooting depth, effectively controlling the rate of root growth. Root parameters for all crops were initialized with R D = 1.0 m, pz = 1, zv = 0 m, and t m = 60 d. With no shallow restrictive layer present at the study site, we assumed the maximum root density was near the surface and therefore zv was fixed at 0, whereas the other parameters were later optimized. Root water uptake in HYDRUS is distributed along the modeled root density distribution, beta, whereas the maximum rate of root water uptake is constrained by the Feddes model (Feddes et al., 2001), described in detail in Kuhl et al. (2018). The parameters of the Feddes model were held fixed at the values provided in the HYDRUS database for maize; the ideal pressure-head range to meet the potential transpiration rate was set between −325 and −600 cm of water, and root water uptake was set to cease below −8,600 cm or above −15 cm of water. Previous results (Kuhl et al., 2018) demonstrated low sensitivity to the lower pressure head limit at which transpiration can be maintained, which could be easily incorporated into the optimization if so desired.
To model major ion chemistry, ionic concentrations of the seven major ions, Ca 2+ , Mg 2+ , Na + , K + , HCO 3 -, SO 4 2-, and Clof the incoming precipitation was specified with the summer seasonal average values (6.0 × 10 −3 , 1.5 × 10 −3 , 0.9 × 10 −3 , 0.4 × 10 −3 , 0.0, 6.5 × 10 −3 , 3.4 × 10 −3 mmol L −1 , respectively) recorded by the National Atmospheric Deposition Program/National Trends Network station MI26 located at the study site. Calcite (CaHCO 3 + ) and dolomite [CaMg(CO 3 ) 2 ] below 2 m were specified as a precipitate in the unleached carbonate zone below 2 m in concentrations of 500 mmol kg −1 each. The coefficient of molecular diffusion in free water and the longitudinal dispersion rates were set at 0.08 cm 2 h −1 and 0 cm −1 , respectively, based on values from Example 2 in the Unsatchem manual (Šimůnek et al., 2012). The soil cation exchange capacity was set to zero based on field data from adjacent sites that indicated the effect on the pore water chemistry due to this process was minor relative to the that of carbonate dissolution . The electrical conductivity of the porewater, σ W (S m −1 ), was calculated by summing the seven major ions, output at each node z (mmol L −1 ), using the ion-specific coefficients from Method 3 proposed by McNeal et al. (1970) for mixed salt solutions: where a n is a linear coefficient and b n is a linear offset for each ion, I n (mmol L −1 ), scaled by the total cation concentration [ ∑ cat =1 ( )]. Separately, we used the PHREEQC software to determine that ion pairing was not an important influence on solution resistivity in the relatively dilute soil porewater solutions above 2 m. An additional 3-yr ramp up period beginning on 1 Nov. 2005 (repeating climate and LAI conditions from 2009 to 2011, data were unavailable prior to 2009) was used for the Unsatchem module runs to stabilize the pore water conductivity, which was very sensitive to the initial ion composition specified as that of the precipitation. Because the root growth and Unsatchem modules could not run in concert, the final root distribution from the root growthenabled model was used with the Unsatchem module. Given the small magnitude of the changes in the modeled pore water conductivity during the growing season in the root zone, we assumed that differences in root water uptake between a static HYDRUS root profile and the dynamic root growth profile had a negligible effect on the ion content of the remaining soil water.

Subsurface electrical resistivity
Subsurface ER is the combined resistivity of the soil-water matrix and is here calculated in a two-step process. First, ER is calculated at a 25˚C reference temperature, then adjusted to in situ temperature. The ER of the soil water matrix, ρ 25 (Ω m), at each node, z, and each survey time, t, was calculated using a modified version of Archie's Law (Archie, 1942) developed by Waxman and Smits (1968): where ϕ [-] is the static porosity, σ S (S m −1 ) is grain surface conductivity, θ (cm 3 cm −3 ) is dynamic modeled water content, σ W (S m −1 ) is modeled conductivity of the pore water, and m and n are dimensionless cementation and saturation exponents. Initial estimates of m and σ S were made by assuming n = 1 and fitting Equation 9 to the soil box-derived petrophysical curve given the porosity and modeled pore water content for each soil layer. These parameter estimates were later updated via optimization with the ER data. Resistivity at in situ soil temperatures ρ T (Ω m) was then calculated using modeled hourly soil temperature, T (˚C), at each depth, z, and time, t. We applied a temperature correction to ρ 25 proposed by Hayley et al. (2007) , valid for temperature ranges between 0 and 25˚C:

Forward geophysical model
The forward geophysical model was built using the boundless electrical resistivity tomography (BERT) model for Python (Rücker et al., 2006). The 2D model domain was 20 m wide by 20 m deep with a fine-resolution triangular mesh (characteristic length = 0.01 m) in the upper 4 m, a coarser mesh (element area = 2 m) to the water table at 17 m, and a very coarse mesh (element area = 5 m) between 17 and 20 m (Supplemental Figure S2). Additionally, the mesh was refined around the location of each buried electrode. The model domain was extended 4 m to the west and east of the 12-m-long ER array, and 18 m deeper than the maximum effective depth to avoid boundary edge effects. The sur-veyed electrode locations, surface topography, and currentpotential electrode pairs for each plot were input to BERT to replicate the dipole-dipole ER surveys (Supplemental Figure S2). For every ER survey date for each plot, the simulated 1D vertical temperature-corrected ER distribution was imported to BERT assuming lateral homogeneity in both soil and root distribution to create a pseudo-2D resistivity distribution.

Optimization
Finally, the coupled hydrogeophysical inversion was accomplished via automated parameter optimization to minimize the root mean sum of squared residuals relative to the TDR and ER data by adjusting the crop and root parameters in the hydrologic model, as well as the petrophysical parameters relating the water content to resistivity. For each experimental plot, our two-step optimization procedure first optimizes hydrologic, root growth, and petrophysical parameters for a single 1D hydrologic model, which is then used to compute laterally-homogeneous ER for the 2D ER model. Then, as a second optimization step, multiple 1D hydrologic and petrophysical models are built, and the soil layer transition depth is optimized for each model. These models each apply to a 1.2-m-wide interval along the ER transect.
Optimization of the first step uses a dual-component objective function, Φ, a weighted root mean square of both resistivity and soil water content residuals of the form: where i is the index of each observation data point in space and time, n is the total number of measurements, θ is the water content, R is the apparent resistivity, and σ is the standard deviation of the dataset. In the first step, Φ was minimized by updating the 16 parameters of interest using a bounded shuffled complex evolution global optimization algorithm (Duan et al., 1992). These parameters included: (a) the crop parameter K C , and root parameters pz, R D , t m for each of the three growing seasons; (b) the petrophysical parameters m for the loam, sandy loam, and loamy sand layers, and (c) σ S for the loam layer only, assuming grain surface conductivity is negligible in the coarser grained textures (Friedman, 2005). The petrophysical parameter n was held constant at 1 for simplification; since this is a lower value than found in the literature, this may produce compensation effects in the esti-mation of the other petrophysical parameters. We elected to not incorporate the soil hydraulic parameters into the optimization routine for several reasons, primarily that the initial texture-based set of parameters described the pregrowing season soil water dynamics quite well, and this greatly reduced the number of unknowns. As it were, global optimization was somewhat resource intensive; optimization of 16 parameters in the coupled model required ∼48 h on a 6-yr old 12-core dual-processor server with 128 GB of RAM. Although not all soils may be so well described by database parameters, adding those with the greatest sensitivity to the optimization would not be prohibitive this approach, as demonstrated in Kuhl et al. (2018). The lower and upper bounds for each parameter were informed by physical constraints and observations from the literature. Parameter constraints and initial estimates are presented in Table 2.
For Equation 11, we first averaged the ER data for each unique effective depth to smooth the along-transect variability in apparent resistivity. This was done for both the modeled and observed resistivity. Soil water content values were taken from the 1D hydrologic model. The RMSEs of the two datasets were weighted by their standard deviations. For soil water content, we used the mean-differenced water content (e.g., θ iobs − θ obs ) to minimize the effect of interlayer variability in soil properties that our simplified layer assumptions did not capture. Parameter estimates were determined by the mean value of the lowest 10th percentile of Φ values from the global optimization. The 95% confidence intervals were calculated as CI 95 = 10 ± 1.96σ 10 √ 10 (12) where 10 , is the mean, and σ 10 , is the standard deviation of each parameter distribution from n 10 iterations within the lowest 10th percentile of Φ values from the global optimization. For our second optimization step, the hydrologic, root growth, and petrophysical parameters were fixed from the first optimization, and only the depth of the transition between the sandy loam and loamy sand soil layers along the transect was estimated. This transition has a strong effect on the apparent resistivity, allowing for lateral sensitivity to the soil properties. Each of the two experimental plots was divided into 10 1.2-m-wide zones, and the depth of the transition between these two layers in each zone was optimized with the full (nonaveraged) 2D ER dataset in a global optimization by minimizing only the resistivity component of Φ. The TDR data were sensitive to any lateral heterogeneity and therefore were not used in the objective function for this portion of the study. Unique root growth parameters for each zone were not estimated with the 2D data.

Vadose Zone Journal
T A B L E 2 Optimization results of each parameter a with initial and final estimates, with 95% confidence intervals, for each plot, as well as the lower and upper parameter bounds. Note that in the miscanthus plot, the m parameter of the loamy sand layer was not optimized due the absence of a substantial loamy sand layer in that plot

F I G U R E 3
Parameter distributions for the 90th percentile (gray) and best 10th percentile (blue) iterations of the global optimization for the miscanthus plot parameters: the crop coefficient K C , root parameter pz, maximum rooting depth R D , and time to maturity t m for each of the three growing seasons, the petrophysical parameter m for the sandy loam and sand soil layers, and the grain surface conductivity σ S for the sandy loam layer only. Similar results from the annual rotation plot are provided in the Supplemental Figure S3. Black dashes represent the 95% confidence intervals, whereas the mean values are reported at the top of each distribution. Greater parameter sensitivity is indicated by a narrower distribution between dashed lines

Parameter estimate distribution
A summary of the results of the parameter estimation via global optimization is presented in Table 2. For both the annual rotation and miscanthus plots, distinct distributions for the best and worst iterations are clearly observed for surface electrical conductivity (σ S ) and the m parameter of the sandy loam soil layer (Figure 3). These parameters have the greatest influence on the objective function Φ, with a clear bimodal distribution in the worst performing iterations distinct from the normal distribution for the highest performing iterations. Estimates of the m parameter varied greatly and were notably estimated to be negative (unrealistic) for the first two soil lay-ers, and near the upper bound for the third and most resistive sand layer. The decision to hold the parameter n fixed at a low value of 1 may have enforced some error into these estimates; however, increasing n would have resulted in an even lower estimate of m due to their direct relationship. For the remaining parameters, the worst iterations are distributed across the full parameter search space, indicating no single parameter estimate is significantly detrimental to the objective function. Among the best iterations, a normal distribution appears for most parameters, indicating convergence towards that value as the best estimate to explain the data. The TDR component of the objective function ranged between 0.5 and 1.1, whereas the ER component ranged between 0.2 and 2.2. There was an overall weak positive correlation between the two components (R = .29 and .35 for the annual rotation and miscanthus plots, respectively) due to the strong influence of the petrophysics on the ER but not the TDR measurements.
Amongst the root growth and crop parameters, the crop coefficient K C , which controls the magnitude of the potential ET, was the most sensitive, with a narrower distribution and confidence interval than the other parameters. The days to maturity parameter t m , which controls how quickly the roots reach R D , was the least sensitive for all years in both plots, with wide confidence intervals and estimates that span the entire search space. The maximum root depth R D was estimated to be greater in the miscanthus than the annual crop in 2009 (2.2 vs. 2.0 m), and then shallower in subsequent years (1.3 vs. 1.7 m) ( Table 2). Notably, 2009, which was the first year after establishment for the miscanthus crop and the year that experienced the most drought-like conditions, shows significant differences in the distribution of the maximum rooting depth, R D , the shape parameter, pz, as well as in K C , compared with 2010 and 2011. Note that increasing R D and decreasing pz together as seen in the 2009 miscanthus estimates has an additive effect, making the ratio of pz/R D smaller, and shifting root fraction downward relative to the other years (Equation 2). In the annual rotation plot by contrast, the covariation in pz and R D results in almost no change from year to year in the ratio of pz/R D . Crop LAI for 2009 and 2011 were based on the most complete year of measurements (2010), and therefore the estimated maximum LAI in 2009 cannot be validated with measurements. Given the importance of LAI in the calculation of potential ET as well as the portioning of E and T (Equations 2-5), the low miscanthus K C parameter estimate for the 2009 growing season relative to the other years could be the model adapting to an overestimation of LAI (and thus potential T) for that period. It also could be reflective of the juvenile nature of the miscanthus crop in its first year post-establishment. We surmise that root depth may be strongly influenced by water content availability regardless of plant age, and that the roots may be deeper, albeit relatively less dense in 2009 than in subsequent years. A lower K C estimate in 2009 may also reflect the physical limitations of reduced root biomass for meeting ET demand. This large effect of LAI on actual ET is most notable in 2009 where canola grew quickly and was harvested early in the season, as well as in the 2011 growing season where the planting of soybean occurred much later than the emergence of the established miscanthus.

Comparison of ET and water content
Modeled outputs of cumulative ET, root depth through time, and water content through time for both cropping systems are shown in Figure 4. Modeled cumulative ET for each growing season (1 May-1 October) in the crop rotation treatment was 304, 395, and 324 mm for canola, maize, and soybean, respectively, compared with ET estimates measured from soil water content loss of 310, 469, and 497 mm for the same plots (Hamilton et al., 2015). Modeled ET for the perennial miscanthus was 351, 430, and 410 mm in 2009, 2010, and 2011, vs. measurements of 516 and 450 mm in 2010 and 2011 (data were not available in 2009) (Figure 4a). Higher modeled ET in the miscanthus plot relative to the annual plot was driven by higher LAI and crop coefficient values and a prolonged season after the annuals were harvested, supporting some previous research findings that miscanthus ET rates may exceed those of annual crops, if soil water is available throughout the growing season. The higher ET modeled in 2010 vs. 2011 in the miscanthus plot was likely due to greater water availability and was also observed by Hamilton et al. (2015). However, our model output for the 2010 growing season was ∼15% lower than the TDR-based estimates from Hamilton et al. (2015) for both the miscanthus and an adjacent maize only treatment plot. Notably, the relative difference between the miscanthus and maize ET was roughly similar in both approaches (46 mm measured difference vs. 35 mm from the inversion). The underestimation of ET ranging between 2 and 34% in our method relative to that used in Hamilton et al. (2015) could indicate some bias in the assumptions for the different models used in both approaches. Modeled water content was not visually much different between the two crops due to the similarities in the timing of the growing season, precipitation rates, and soil properties that largely constrain the water balance. Modeled growing season ET differed by between 8 and 20% between the two plots over the 3 yr.
Using ROSETTA parameters to model soil water retention, the absolute value and range of the modeled water content and the temporal dynamics show good agreement with the observations in the sandy loam (0.2-m depth) and sand (1.25-m depth) soils for both crops (Figure 4c-e). The largest nonwinter discrepancy between the model and the deeper soil moisture measurements occurs in the miscanthus crop in the fall seasons of 2009 and 2011 where the first large pulse of percolating water after the growing season is underpredicted (Figure 4e). In the shallower water content measurements, most precipitation pulses are well modeled with two exceptions: a large wetting event in August 2011 that coincides with the driest conditions and causes a large pulse of infiltrating water in the observations for both plots but not the model; and notably in the fall of 2010, where observed miscanthus water content does not increase with precipitation events as seen in the annual rotation nor in the model of either crop. This may be due to substantial crop residue at the surface that limits infiltration. Aside from this, agreement in the other four depths is also good in both field plots, with low RMSE (0.027 for miscanthus and 0.022 cm 3 cm −3 for the annual rotation) across 3,510 and 3,660 total observations for each respective crop spanning the three model years (Figure 4d

Comparison of modeled and observed root depth
Model parameterization of root growth and maximum depth through time (Figure 5b) yielded reasonable estimates for both crops, with root fraction decreasing to <1% of the maximum root density at depths of 1.32, 1.20, and 1.24 m for canola, maize, and soybean respectively, and 1.64, 0.80, and 0.72 for miscanthus for each model year. Estimates for the annual rotation fall within the ranges of temperate agricultural crops summarized in Fan et al. (2016) . The estimated root depths for miscanthus agree with the shallow root distribution for rainfed miscanthus found by Mann et al. (2013); however, this is shallower than the roots observed to depths of 2.7 m by Ferchaud et al. (2015) and 1.8 m by (Neukirchen et al., 1999). Note that the root depth shown in Figure 5b is much shallower than the parameterized R D values shown in Table 2 and Supplemental Figures 3 and S3, which represent the depth where root fraction becomes 0. The exponential decay of root density (Equation 2) produces a long tail of nominal root fraction values which we chose to truncate at 1%. We hypothesize that the deeper root depth for both crops in 2009 is likely a response to the drought conditions, which is supported by the limited measurements of root biomass collected at the end of the growing season using soil cores at the center and adjacent to the plant, which show a deeper distribution of normalized root mass in 2009 than in either 2010 or 2011 (Figure 5b). The model generally overestimates the root fraction (normalized by the maximum model and measurement values) with depth relative to the data ( Figure 5b); however, this could be explained by the decay of deeper fine roots between the growing season and the sampling period after senescence in December. No root samples were collected from below 1 m, so we cannot compare maximum rooting depth; however, given the low root density in the 0.5-to 1.0-m core length (Figure 5a), we infer that relatively few roots grew to depths below 1 m at this location. Deeper root distribution does not directly lead to greater ET across crops in our simulations, as maximum root depths were found to be greater for the maize and soybean than miscanthus in 2010 and 2011, respectively. Although in situ root distribution data for validation were only available for miscanthus, the modeled root distribution estimates generally agree with F I G U R E 5 Depth profiles of (a) average log observed root mass per volume with 95% confidence interval (CI) error bars (n = 5) from four depth cores (shaded region indicates length of each core) sampled at the end of each growing season for miscanthus; (b) observed root mass per volume normalized by the maximum measurement for each of the three sample dates (open circles) and modeled normalized root fraction distribution (dashed lines) at the end of the growing season for each year. Note that the 2009 sample was not collected until the following spring, and the downward shift in root distribution could be attributed to root decomposition in the upper 0.1 m after the end of the growing season. 2009 was also the first year after establishment from rhizomes, which could be expected to have less root mass than subsequent years the data, where root mass was found to be concentrated in the upper 0.5 m.

Hydrogeophysical model outputs
In addition to model performance through time, we show for the 2010 growing season the 1D hydrogeophysical model outputs and observations, where available, of root distribution, water content, temperature, and porewater electrical conductivity (inverse of resistivity), all of which are used to calculate the temperature-corrected ER with depth via Equations 9 and 10 for each ER survey date ( Figure 6). Overall, the model reasonably simulates the dominant trend of decreasing water content and increasing temperature throughout the growing season and the resultant resistivity changes. Modeled porewater ion content, and hence conductivity, varies with depth and time but over a range of only about ±5% (Figure 6d), resulting in smaller percent effects on the bulk soil resistivity (Equation 9). Temporal changes in soil water content and temperature, which tend to be correlated, are the dominant factors driving observed resistivity changes, with the clearest increase in resistivity occurring between the wetter early season surveys (blues and green) and the drier growing season surveys (yellow and pink) (Figure 6e). Amongst the early surveys, the most resistive in both the model and the observations is the earliest survey on 3 May 2010, due to the much colder temperatures (Figure 6c), despite a slightly higher pore water conductivity (Figure 6d), which has the opposite effect. The largest error in the modeled resistivity occurs across the four largest electrode spacing geometries (e.g., those with the deepest prospecting depth), which could be indicative of limitations of the petrophysical relationship introduced by fixed parameters or the pore water conductivity model. It also could reflect error in the deeper (>1.25 m) water content estimates, which were beyond the available depth of the TDR data.

Modeled and observed change in resistivity
To summarize how well the model predicts the changes in resistivity due to hydrogeophysical processes across the two plots and three growing seasons, we calculated the change in apparent resistivity for each unique geometry and compared the model to the observations (Figure 7). Generally, the modeled magnitude and direction of change agree approximately with the observations for both plots, particularly later in the growing season of each year during extended drying periods, which the model predicts quite well (Figure 7). We also note that the change in resistivity is generally similar across the two crops, although there are notable exceptions. In 2009, there is a stark difference between the 1 June and 29 June 2009 surveys, where the observations in the annual rotation plot data show an increase in resistivity while miscanthus shows F I G U R E 6 One-dimensional (1D) outputs of the 2010 growing season for the miscanthus plot, showing the model results (solid line) and observations (circular markers) at monthly intervals from May to September, for (a) normalized root density distribution, (b) volumetric water content (VWC), (c) temperature, (d) porewater electrical conductivity (EC), (e) log resistivity, and (f) apparent resistivity in the upper 2 m of the model domain. The depth of the apparent resistivity curve is the effective depth (a function of electrode spacing) and is only an approximation. The strong contrast in volumetric water content and surface electrical conductivity due to the presence of clays between the sandy loam and loamy sand soil interface at ∼0.5 m drives the shape of the modeled resistivity and apparent resistivity curves F I G U R E 7 Summary of modeled (Mod) vs. observed (Obs) changes in the one-dimensional (1D) apparent resistivity between consecutive surveys in each model year for the annual rotation (open circles) and the miscanthus (stars) plots. Positive delta values indicate an increase in resistivity with time due to drier conditions during that period, whereas negative values indicate a decrease in resistivity with time due to wetter conditions. The amount of change scales with the absolute value of apparent resistivity, which increases with depth a decrease (Figure 7a). This discrepancy is not seen in the model, which predicts a decrease in resistivity for both plots due to a large precipitation event on 19 June 2009. Similarly, the change between the 28 June and 26 July 2010 surveys is undermodeled, with data for both plots showing an increase in resistivity while the model predicts a small decrease (Figure 7b). This is likely due to a large precipitation event occurring on 22 July 2010 that caused a greater pulse of infiltrating water in the hydrological model than is observed in the data.

3.2.3
Estimate of sandy loam interface Using the estimated 1D hydrogeophysical models for each crop to estimate the spatial heterogeneity in the depth of the bottom of the sandy loam soil layer resulted in good agreement in the 2D ER data (Figure 8). The trend in depths of this layer interface roughly follows the surface topography (Supplemental Figure S2) in both plots, indicating a relatively flat laying layer that is deepest at the highest elevation to the west and becomes shallower with decreasing surface elevation to the east, with a slight rise in both plots west of the center. Depth estimates ranged from 0.50 m in the first zone to 0.38 in the 10th zone with a mean of 0.48 m for the annual rotation plot and from 0.52 m in the fourth zone to 0.38 in the eighth zone with a mean of 0.45 for the miscanthus plot (Table 3). This variability in the depth of this interface at the plot scale seems reasonable given the similar range of interface depths observed in core samples at each of the treatment plots, including those not included in this study. At these 10 sample points, the depth of the interface ranged from 0.3 to 0.7 m. For both crops, the modeled ER was generally overestimated closer to the surface where resistivity was low, and underestimated at measurements with the larger electrode spacings, where resistivity was high (Figure 8e, f), suggesting that our simplified assumption of three sets of petrophysical parameters for the three primary soil textures may not be sufficient to capture the vertical heterogeneity. It should also be noted that the optimized petrophysical curve for the sandy loam soil was considerably lower than the soil box resistivity model; modeled wilting point resistivity values did not exceed 100 Ω m vs. the 255 Ω m observed in the soil box experiments (Table 1). Conversely, the modeled resistivity ranges for the sand layer were roughly three times greater than the soil box resistivity (Supplemental Figure S1). Above the carbonate zone that begins at 2 m depth, the only ionic input to the pore water conductivity model came from dilute precipitation, whereas other sources of ionic input (particularly nitrate leaching in the fertilized maize and mineral weathering of aluminosilicates) that could modestly increase pore water conductivity were not included. Error in the mod-T A B L E 3 Estimated transition depth to coarser soil texture with lower water holding capacity in each of the 10 zones from 0 to 12 m (west to east; Supplemental Figure S2  F I G U R E 8 Two-dimensional (2D) pseudosections (shown here without topographic correction) of (a) observed and (b) modeled apparent resistivity from the annual rotation plot and (c) observed and (d) modeled apparent resistivity from the miscanthus plot on 24 Aug. 2009 after spatial optimization of the transition depth at the bottom of the sandy loam soil layer. White gaps are data outliers that were filtered and neither modeled nor included in the objective function calculation. Modeled vs. observed ER for all 17 survey days for the (e) annual rotation, and (f) miscanthus crops. Note that overall ER is higher for the miscanthus plot than the annual rotation, due to the absence of an intermediate loamy sand soil and a shallower transition to sand in the miscanthus plot eled pore water conductivity linearly propagates to the resistivity estimates and could significantly alter the estimated petrophysical parameters, although the estimated m parameter can accommodate for some of this error. Additionally, although a greater number of zones may have resulted in better agreement between the model and observed ER, for the purposes of this study, we limited the number of zones to 10 for efficiency. The petrophysical relationship between water content and resistivity remains one of the largest sources of uncertainty in this modeling approach and further constraints on the shallow pore water conductivity dynamics would be advantageous to build confidence in the major ion chemistry modeling work proposed here.

Sensitivity of ET to soil interface depth
To analyze the sensitivity of seasonal ET to realistic variations in depth of the interface between high and low field capacity soils, we compared the cumulative ET in each zone to the 1D model with a fixed interface depth at the 10-zone average of 0.48 and 0.45 m for the two cropping systems, respectively (Supplemental Figure S4). There is a clear correlation between drought conditions and ET discrepancy, with the largest difference in 2009 and the smallest in the beginning of 2010, when there was consistent rainfall throughout the growing season (Supplemental Figure S4). When high ET demand draws down soil water content (i.e., during dry periods), the presence of a deeper sandy loam soil with more plant available water becomes more important.  Figure S4). However, under the drier conditions in 2011, the miscanthus crop shows a much larger decrease in ET than the soybean crop, likely due to the roots being shallower (Figure 4), which limits access to deeper water reserves.

CONCLUSIONS
Here, we demonstrate the use of a broadly applicable novel coupled hydrogeophysical inversion approach to estimate root growth and root water uptake distribution properties, along with heterogeneous soil properties below nonirrigated biofuel crops in a field setting. With monthly ER surveys over three growing seasons and continuous hourly TDR data, we estimated the crop parameters of multiple hydrological models to describe the root water dynamics of a canola-maizesoybean annual crop rotation and a perennial miscanthus grass crop. Using HYDRUS-1D, we built for the first time to our knowledge, a forward process-based model that simulates root growth, water content, pore water conductivity, and temperature to predict changes in ER through time, using a petrophysical relationship with estimated parameters. We found good agreement between modeled and observed apparent resistivity, as well as the change in resistivity through time across all 17 ER surveys in both crop treatments on the sandy loam soils at our study site. Although the simulated water content reasonably matched the TDR data throughout the model period at all six observation depths, during dry periods, the magnitudes of pulses of infiltrating water were poorly estimated by the hydrological model, resulting in greater discrepancies between the ER model and observations near those events. On average, cumulative ET was reasonably estimated from the coupled inversion; discrepancies with data-driven estimates (Hamilton et al., 2015) varied from excellent agreement (2%) to significant underestimation (34%). The perennial cropping system consistently had higher ET than the annual crop rotation, driven by a longer growing season and higher LAI, resulting in a larger fraction of potential ET as potential transpiration, and a greater estimated crop coefficient that contributed to higher potential ET. Reasonable estimates of root depth for all crops were also achieved, with the method predicting deeper roots in 2009 during the most drought-like conditions. The rooting depth did not correlate with the modeled ET, as miscanthus had a shallower root distribution than maize and soybean in the 2010 and 2011 model years. The deepest modeled root distribution was estimated for miscanthus in 2009, which may suggest the grass has a favorable adaptation to drought conditions, with implications for planting on marginal lands. Using the calibrated hydrological models and the 2D ER data to estimate the field-scale heterogeneity in soils, we resolved the lateral variation in thickness of the upper sandy loam soil layer from 0.38 to 0.57 m. We found the relatively high water availability in this uppermost relatively loamy soil layer to have a moderate influence on the cumulative growing season ET, particularly in 2009, when the driest conditions were observed.
The ER data are very sensitive to changes in soil water content, and large transects including three-dimensional regions can be surveyed with minimal soil disturbance. This study demonstrates that coupling ER data with TDR soil water content data and a process-based hydrological model allows interpolation temporally between ER surveys and extrapolation beyond the point scale of soil water content sensors. While refining the estimated petrophysical parameters with the coupled inversion enabled improved fit between the data and the model, this approach yielded smaller than expected values of the cementation (m) exponent of Archie's Law, likely due to compensation for errors in the pore water conductivity model. Improvement of petrophysical estimates could be made in the future with detailed shallow pore water conductivity data to aid in the calibration of dynamic relationship between water content and resistivity. This method worked well at this Michi-gan study site with relatively dry, high-resistivity loamy sand soils and should be applicable to a variety of settings, provided some contrast in spatial and temporal water content variability exists during the period of interest. Although we chose to use database parameters for the soil water retention function, this method is well suited to estimate those parameters in the joint optimization of geophysical and hydrological data as demonstrated in Kuhl et al. (2018). The presented approach holds promise as a tool to better understand root water dynamics of crops and other vegetation, improving our ability to make inferences about the implications of land use change on regional water balances and informing cropping decisions including the transition from annual crops to large perennial grasses. As we look towards a clean energy future, understanding and minimizing the environmental impacts of large-scale conversion of land to biofuel crops will remain a critical research question.

A C K N O W L E D G M E N T S
This material is based upon work supported in part by the Great Lakes Bioenergy Research Center, the USDOE, Office of Science, Office of Biological and Environmental Research under Awards no. DE-SC0018409 and DE-FC02-07ER64494, the National Science Foundation Long-term Ecological Research Program (DEB 1832042) at the Kellogg Biological Station, by Michigan State University AgBioResearch, and by the USDA's National Institute of Food and Agriculture Award 2015-68007-23133. This work was made possible by the efforts and support of many. In particular, we would like to thank Kaya Diker for his extensive work on the electrical resistivity array installation and data collection; Kevin Kahmark and Mir Zaman Hussain for assistance with collecting the soil cores; Randall Schaetzl for the grain size analysis and complimentary use of his soil laboratory; and Amanda Liddle for her effort preparing soil samples and conducting the soil box resistivity tests.