Rothamsted Repository Download

Geophysical methods, such as electromagnetic induction (EMI), can be effective for monitoring changes in soil moisture at the field scale, particularly in agricultural appli-cations. The electrical conductivity ( σ ) inferred from EMI needs to be converted to soil moisture content ( θ ) using an appropriate relationship. Typically, a single global relationship is applied to an entire agricultural field; however, soil heterogeneity at the field scale may limit the effectiveness of such an approach. One application area that may suffer from such an effect is crop phenotyping. Selecting crop varieties based on their root traits is important for crop breeding and maximizing yield. Hence, high-throughput tools for phenotyping the root system architecture and activity at the field scale are needed. Water uptake is a major root activity and, under appropriate con-ditions, can be approximated by measuring changes in soil moisture from time-lapse geophysical surveys. We examine here the effect of heterogeneity in the θ – σ relationship using a crop phenotyping study for illustration. In this study, the θ – σ relationship was found to vary substantially across a field site. To account for this, we propose a range of local (plot specific) θ – σ models. We show that the large number of parameters required for these models can be estimated from baseline σ and θ measurements. Finally, we compare the use of global (field scale) and local (plot scale) models with respect to ranking varieties based on the estimated soil moisture content change.


INTRODUCTION
Over the past two decades, there has been a growth in the use of geophysical methods in agriculture (Allred, Daniels, & Ehsani, 2008). This has been driven, in part, by the need to assess variation in soil properties in a noninvasive manner over relatively large scales. Geophysical methods in such a context are a subset of proximal soil sensing approaches (Viscarra Rossel, Adamchuk, Sudduth, McKenzie, & Lobsey, 2011). Measurements of properties, such as electrical conductivity, are typically treated as a proxy for a soil property or state of interest (e.g., soil texture, bulk density, or soil moisture content). Such methods may also be used in a time-lapse manner to examine changes in soil properties or states (e.g., changes in texture or soil density due to land management practices). Typically, maps of a geophysical property are presented in a qualitative manner. Although this can be effective in some cases, the ability to estimate quantitatively the property, or state, of interest offers greater scope for a wider range of agricultural applications. To achieve such quantification, the relationship between the geophysical proxy and the soil property or state is required. Such relationships may be spatially variable, particularly over field scales typical in agricultural studies. Here, we assess such heterogeneity in a wheat (Triticum aestivum L.) phenotyping study and propose practical methods to account for such variability.

Field-scale phenotyping bottleneck
Wheat is one of the main staple crops in the world. It has been bred over centuries for specific traits, most of which are aboveground characteristics. Given uncertain future climatic conditions, there are demands for more resilient breeds.
A key component of such resilience lies in the root system of the crop. Deeper root systems are correlated with higher yield and higher resistance to drought (Wasson et al., 2012). Usually the root system of a crop is assessed in the laboratory or in the greenhouse. However, field studies of the root system are essential to understand more about how each variety adapts to its environment. The typical approach of assessing the root system of a crop in the field is by direct sampling (Wasson et al., 2014). Such methods are destructive, labor-intensive, and expensive in a conventional breeding program with a large number of breeding lines. An alternative, less invasive, and quicker approach is to consider the root activity rather than the quantity of roots. Such methods rely on observing changes in soil moisture to infer root activity (Beff, Günther, Vandoorne, Couvreur, & Javaux, 2013;Garré et al., 2013;Michot et al., 2003;Srayeddin & Doussan, 2009). Different methods to measure efficiently this change in soil moisture were explored by Whalley et al. (2017) for different wheat genotypes. Among them, geophysical methods, such as electrical resistivity tomography (ERT) and electromagnetic induction (EMI) appear promising as a means of measuring a proxy to observe the dynamics of soil moisture of the subsurface . Shanahan, Binley, Whalley, and Watts (2015) illustrate the use of EMI for differentiating soil drying from different wheat genotypes in a phenotyping context. In their study, the relationship between the observed proxy (soil apparent electrical conductivity) and soil moisture content was assumed to be homogeneous across the study site. Huang, Purushothaman, McBratney, and Bramley (2018) also use EMI as a proxy for plot-scale crop water of different chickpea (Cicer arietinum L.) genotypes. Other examples of the use of EMI in crop-related studies include Cassiani et al. (2012), von Hebel et al. (2014, and Moghadas, Jadoon, and McCabe (2017).

Electromagnetic induction
The EMI method measures the soil apparent electrical conductivity (σ a ) in a noncontact or invasive manner. A standard Core Ideas • Field-scale relationships between θ and electrical conductivity can be inappropriate. • Pedophysical parameters can, in some cases, be approximated using baseline data. • The method is illustrated for application of EMI mapping for phenotyping wheat crops.
EMI device is composed of a transmitter (Tx) coil and at least one receiver (Rx) coil. The transmitter coil generates a transient electromagnetic field. This primary field induces eddy currents in the ground; the magnitude of eddy currents generated is a function of the soil electrical conductivity, σ. The eddy currents then induce a secondary electromagnetic field. Both primary and secondary electromagnetic fields are measured by the receiver coils. The out-of-phase component of their complex ratio is used to compute the apparent electrical conductivity (σ a ) of the subsurface. Electromagnetic induction measurements can be made in vertical and horizontal coplanar orientations, with different depth-sensitivity functions. Several current instruments, such as the one used in this study (Mini-Explorer from GF-Instruments), have multiple receiver coils. The relationship between depth-specific σ and measured σ a , for a given coil orientation, and the distance between the Tx and Rx, can be described using a simple function-the "cumulative sensitivity function" (McNeill, 1980). A more accurate, but more complex, method based on Maxwell's equations (Andrade, Fischer, & Valenta, 2016;von Hebel et al., 2014) can also be used to describe such a relationship. Using measurements made on a multi-coil device, depthspecific σ can be determined from inverse modelling of the σ-σ a relationship. The inversion process seeks the best distribution of depth-specific σ that is consistent with all observed σ a values for different coil spacings and orientations. A prerequisite, considered by some authors, for inversion is that the apparent values given by the different EMI configurations need to be calibrated with results from an ERT survey (Lavoué et al., 2010). More details about EMI inversion can be found in von Hebel et al. (2014).
Electromagnetic induction measurements have been extensively used to map field heterogeneities and produce detailed soil maps for the definition of management zones in precision agriculture (Brevik, Fenton, & Lazari, 2006;Corwin & Lesch, 2003;King et al., 2005). More recently, multi-coil EMI instruments have provided greater depth-specific information in agricultural studies, allowing assessments of depth-specific σ and its link to aboveground crop performance indicators (Brogi et al., 2019;von Hebel et al., 2018).

Soil moisture content-electrical conductivity relationships
The soil electrical conductivity is controlled by a number of properties (soil texture, organic matter content) and states (soil temperature, pore water electrical conductivity, bulk density, and soil moisture content). The soil structural state and its properties control σ through pore connectivity and porosity. Such properties are also inherently linked to soil moisture content (e.g., determining residual moisture content), which has a major effect on soil σ. Temperature effects can be accounted for given local vertical soil temperature profiles, which we assume to not vary spatially inside the same field, although effects of daily or seasonal variation in temperature may need to be accounted for. The electrical conductivity of the pore water also contributes to the soil σ. In temperate climates, the variation of the pore water electrical conductivity should be minimal in rain-fed settings. However, this has a greater impact in irrigated conditions, as the irrigated water (e.g., groundwater sourced) is likely to have a different ionic composition and temperature than the pore water in the surface layers of soil. In semiarid environments, pore water conductivity effects may be significant due to enhanced salinity arising from high evaporative fluxes (Corwin & Lesch, 2005). Note that even in rain-fed environments, increase in pore-water electrical conductivity can occur due to fertilizer application.
Archie's law (Archie, 1942), developed for oil reservoir investigations, is a commonly used empirically derived model that relates the soil condition to the bulk σ. Waxman and Smits (1968) extended Archie's law by accounting for the effect of clay minerals (forming surface electrical conductivity). Several other approaches have been developed specifically for soils (Rhoades, Raats, & Prather, 1976). Laloy, Javaux, Vanclooster, Roisin, and Bielders (2011) compared a range of models for soil electrical conductivity, adopting the term "pedo-electrical" model to differentiate this from classical petrophysical approaches.
Following Laloy et al. (2011), the relationship between σ and soil moisture content (θ) can be expressed as where a, b, and n are empirical parameters that depend on soil properties. Following Garré, Javaux, Vanderborght, Pagès, and Vereecken (2011), a is influenced by the pore water conductivity, soil texture, and porosity; b is influenced by the soil surface conductivity; and n is controlled by the soil texture. When the exponent n is close to 1, Equation 1 can be approximated by a linear relationship. The parameters of Equation 1 may be obtained from laboratory measurements on field samples (Shanahan et al., 2015) or directly in the field, for example using a trench and soil moisture sensors (Beff et al., 2013;Garré et al., 2013;Michot et al., 2003). Both methods provide information on a relatively small volume that might not be representative of the entire field. Indeed, from field-scale observations, the different soil textural properties also affect the θ-σ relationships, either when using σ a (Stanley, Lamb, Falzon, & Schneider, 2014) or with depth-specific σ (Jayawickreme, Van Dam, & Hyndman, 2010). Equation 1 is usually appropriate when the soil moisture change is large and the soil heterogeneity is small. However, if significant soil heterogeneity exists, the variation in the parameters in Equation 1 may need to be accounted for. This effect may be particularly important in phenotyping studies (the determination of specific traits of crop varieties), since the differences in soil moisture change between crop lines (varieties) may be smaller compared with other studies where different species are used. Whether depth-specific or apparent values (like in this study) are considered, estimates of small changes in soil moisture are likely to be affected by heterogeneity in the θ-σ relationship (Equation 1).
Furthermore, in a phenotyping context, a better prediction of the soil moisture or change in soil moisture from EMI is important, as it can help to make the variety ranking similar to the one obtained with direct soil moisture observations. Of course, if direct soil moisture data are available, there is little value in additional geophysical proxy measurements. However, in this study, the direct measurements allow us to determine what the maximum achievable information on soil moisture content obtainable from EMI measurements might be.
Therefore, this study aims (a) to quantify the spatial heterogeneity of θ-σ relationships at the field scale, (b) to determine its impact on the phenotype ranking of wheat lines, and (c) to explore approaches to account for such effects using simplified but practical approaches. The investigation uses a dataset of σ and θ measurements collected during a winter wheat field experiment.

Field layout
Measurements were made during the 2016-2017 growing season at the Warren Field experimental farm (Woburn, UK; 52 • 01 ′ 06.5" N, 0 • 35 ′ 29.0 ″ W) operated by Rothamsted Research. The soil at the site is classified as a sandy clay loam (Distric Cambisol with 54% sand, 20% silt. and 26% clay, more details in Shanahan et al., 2015). The field was sown with winter wheat at the end of 2016 and harvested in August 2017 (Bai et al., 2019). In the experiment, 71 lines of wheat and one fallow treatment (all with three replicates) were randomly distributed in three blocks. An aerial photograph showing the field experiment and the 216 plots is shown in Figure 1. Out of the 216 plots (each 9 by 1.8 m), 12 plots (four F I G U R E 1 Aerial picture of the field showing the 216 plots (each 9 × 1.8 m) sown with winter wheat in 2016. Plots marked in red are equipped with electrical resistivity tomography (ERT) arrays varieties) were equipped with a 24-ERT array (0.25-m spacing) placed along the middle of each plot. The ERT data were used to calibrate EMI measurements following Lavoué et al. (2010). All plots were equipped with a 1.5-m-long neutron probe access tube positioned 1 m from the edge of the plot. The ratio counts from the neutron probe were converted to soil moisture content using a field calibration (±0.01 cm 3 cm −3 ). In the field, temperature sensors recorded soil temperature at 0.1-, 0.2-, 0.3-, 0.4-, 0.6-, and 1-m depths. They were used to correct the electrical conductivity from the ERT and EMI using the ratio model (Ma, McBratney, Whelan, Minasny, & Short, 2011) with a 2% increase per degree Celsius.

Field measurements
Three sets of EMI measurements were collected on each plot with a Mini-Explorer instrument (GF Instruments) according to the guidelines provided in Shanahan et al. (2015). They were then averaged to obtain a mean for each plot. Surveys were conducted on dates: 8 Oct. 2016, 2 Mar. 2017, 16 Mar. 2017, 3 Apr. 2017, 27 Apr. 2017 May 2017, and 1 June 2017. Data from some plots were discarded because of two high voltage cables buried under the field. The filtering used the standard deviation of the three sets of EMI data for each plot.
The Mini-Explorer contains three receiver coils with separations 0.32, 0.71, and 1.18 m from the transmitter coil. Measurements in the two modes (horizontal coplanar mode [HCP] and vertical coplanar mode [VCP]) were obtained. Therefore, six measurements of apparent conductivity were made. The normalized sensitivity pattern (McNeill, 1980) of each configuration is shown in Figure 2a (note that in Figure 2a and hereafter, the notation [e.g., HCP0.32] is used to identify coil orientation and spacing [HCP with a 0.32 m coil spacing]). F I G U R E 2 (a) Normalized local sensitivity pattern for the six pairs of coil orientations and coil separations available on the Mini-Explorer instrument. The triangles show the depth above which there is 70% cumulative sensitivity (commonly referred to as the effective depth of investigation). (b) Measured soil moisture content profile by neutron probe. To build the apparent soil moisture content, each depth-specific soil moisture content (θ) measurement is multiplied by the integrated electromagnetic induction (EMI) sensitivity corresponding to its depths (between the gray lines) and then summed (see Section 3.1). HCP, horizontal coplanar mode; VCP, vertical coplanar mode Figure 2b shows example soil moisture data from the neutron probe taken at seven depths. For each depth, the gray lines denote the limits used to compute the local sensitivity weights used in the computation of the apparent soil moisture content (Section 3.1).
The ERT measurements were collected using a 48 Syscal Pro (Iris Instrument)  Whenever possible, ERT and EMI measurements were collected on the same day. Neutron probe datasets were collected as close as possible to the ERT/EMI dataset, either on the same day or before or after an interval of a few days, thus minimizing disturbance from any rainfall events. Note that the neutron probe dataset of mid-May was taken after a large overnight rainfall event. This had an impact on the shallow measurements (0.15and 0.30-m depths) but did not influence the deeper measurements. Note also that N fertilizer was applied just before the measurement at the end of May. However, because of its F I G U R E 3 (a) Rainfall and potential soil moisture deficit (PSMD) with markers corresponding to the collection date of the electrical resistivity tomography (ERT), electromagnetic induction (EMI), and neutron probe (NP) dataset. (b) Evolution of soil apparent electrical conductivity (σ a ) from EMI. (c) Evolution of computed apparent soil moisture content (θ a ). (d) Evolution of the measured soil moisture content from neutron probe for selected depths. Error bars are SEM (sometimes too small to be visible on the graph). Dotted lines are averages of the fallow plots, whereas solid lines are averages of the cropped plots. HCP, horizontal coplanar mode; VCP, vertical coplanar mode application as dry pellets and the lack of large rainfall events, it is unlikely that it had fully dissolved into the soil at the time of the end of May survey. This could have caused a significant increase in the pore water electrical conductivity, and hence in our EMI and ERT measurement, no sharp increase in observed electrical conductivity is apparent. At the end of the field campaign, four different datasets of ERT, EMI, and neutron probe measurements were available to derive pedophysical relationship for each plot. Despite the limited number of time-lapse data collected on the same plot, the larger number of plots screened enables us to capture well the temporal and spatial variability across the field.

Apparent soil moisture content
To allow comparison with observed apparent conductivity measurements and to avoid any inversion artifacts that can arise from EMI inversion, an "apparent" soil moisture was computed based on the weights of the EMI cumulative sensitivity function (Figure 2a) following the approach given by Martini et al. (2017). The θ measurements of a given profile ( Figure 2b) were multiplied by their respective depth-specific normalized local sensitivity and then summed to obtain an apparent soil moisture content (θ a ). The shape of the normal-ized sensitivity function is determined by the same parameters as for the EMI: the coil orientation (HCP or VCP) and the coil spacing (0.32, 0.71, or 1.18 m). Thus, for each pair of coil orientation and coil spacing, a different θ a was obtained, for comparison with the observed σ a from EMI. The apparent soil moisture content θ a is given by where θ i is the measured soil moisture content of layer i, and s i is the sensitivity of the layer i derived by integrating the cumulative sensitivity function between the top and the bottom depths of the layer (Figure 2). Note that the sum of s i for the profile is equal to 1. n is the number of layers.

Time-lapse approach
Time-lapse monitoring of σ allows the removal of stationary effects of the soil (soil organic matter, soil texture) on the θ-σ relationship (Robinson, Abdu, Lebron, & Jones, 2012;Shanahan et al., 2015). This approach relies on the measurements of a baseline (in this case, where no crop effect is present), which is usually made at the beginning of the growth season when the field is at or near field capacity. All subsequent surveys can be compared with this baseline, consequently revealing the main drying pattern mainly driven by root activity. For the experiment presented here, the baseline data were measured on 16 Mar. 2017. There are two ways to compute changes from the baseline conditions: (a) by computing the difference, or (b) by computing the relative change. Assuming a linear relationship between θ and σ (n = 1 in Equation 1), the equations below can be written.
The difference is simply the difference between σ and σ ref : where σ ref and θ ref are the baseline σ and θ, respectively. The relative change is the difference between σ 1 and σ ref normalized by the baseline σ ref (Equation 4). It is given by Computing differences (Equation 3) removes the effect of "offset" b but retains "slope" a, which may vary across the site. In contrast, working with relative change (Equation 4) retains the effects of a and b, unless b is relatively small. In the latter case, Equation 4 can clearly be simplified to link directly the relative change in σ with the relative change in θ as The expressions above were used to explore ways in which the variation of a and b within a site can be accounted for. Figure 4 shows the different relationships between σ a and θ a for three plots with the same variety in the field site. The variation between the three responses (expressed as absolute, difference, or relative change) reveals the effect of spatial variability across the site, highlighting the limitation of adopting a single global relationship. Figure 5 shows the distribution of θ a and σ a in April 2017 and their difference with respect to the baseline in March 2017 (16 Mar. 2017). From Figure 5, it can be seen that the patterns for both absolute and differences are different. This illustrates the effect of different θ-σ relationships observed in Figure 4. Both patterns in σ a and θ a values remain consistent for the different collection dates.

Development of local model
Typically, a few samples from the field are collected to build a global unique relationship between θ and σ. We can express this relationship as σ = g θ + g (6) where the global a g and b g parameters are identical for all the plots. However, for a heterogeneous field, using this global relationship may lead to substantial errors in the estimation of soil moisture content changes. In order to overcome this, we explored local models allowing the assignment of a unique θ-σ relationship for each plot: M1. Linear local model, based on Equation 1 assuming n = 1. This model has two plot-specific parameters: i is the plot number, the slope is a i , and the offset is b i : Figure 6 illustrates, using all measurements, how well the linear global model and linear local model (M1) perform. There is a clear (and expected) improvement of the prediction of soil moisture content with the linear local model. Note that an exponential model (not shown here) following Equation 1 was also fitted and has similar performance to the linear model (R 2 = .37 for the global exponential model; R 2 = .82 for the local exponential model). Consequently, the linear model is adopted hereafter.
As seen in Figure 6, the local linear model outperforms the global linear model but increases the number of parameters needed. More importantly, a full set of monitored soil moisture content values is needed, making the geophysical proxy approach redundant. As a first step to reduce the number of local parameters, we introduce two new models: M2. Multi-offsets model: a linear model where each plot has its own offset b i but share a common slope a g : M3. Multi-slopes model: this model only applies to differences in values and is based on Equation 3, with each plot having its own slope a i . This model has one parameter per plot (slope): Mathematically, the multi-offsets model (M2) produces a set of parallel σ-θ relationships similar to Figure 4a, whereas the multi-slopes model leads to a set of conical Δσ-Δθ relationships similar to Figure 4b. Both use fewer parameters than the local linear model (M1). The rationale for these simpler models is the need to reduce the number of parameters needed and increase our ability to predict them using a set of baseline measurements.

Development of predicted local (plocal) models
All local models (M1-M3) require large amount of information for each plot and have limited practical use in a field F I G U R E 6 Both graphs show the observed apparent soil moisture content (θ a ) vs. the predicted θ a from (a) the global linear model (Equation 6) and (b) the local linear model (Equation 8) phenotyping application. As stated above, if direct measurements of soil water were available in a field experiment, there would be no benefit or value in using alternative geophysical proxy measurements. However, they allow us to determine what the maximum achievable information on soil moisture content obtainable from EMI measurements might be. As a more practical solution, we explore a range of alternative approaches where the local θ-σ relationship is known for a subset of plots and the geophysical data are used to predict those local relationships for the other plots (plocal).

Predictors of the local parameters
The first step in developing predicted local (plocal) models is to identify the best estimates of the local parameters among baseline measurements. Figure 7 shows the relation-ship between the different local parameters from each model (M1-M3) and the baseline σ a and θ a . It can be observed for the linear local model (M1) that the local offsets (b i ) are well related to baseline θ a ref and that the slopes (a i ) are more related to σ a ref .
The multi-offsets (M2) and multi-slopes (M3) models aim to amplify those trends by reducing the number of local parameters. Using multiple local offsets but a global slope (Equation 9), the multi-offsets model (M2) displays a stronger relationship with the baseline θ a ref (R 2 = .86) than the linear model (R 2 = .40). Using multiple local slopes and no offsets (Equation 10), the multi-slopes model (M3) displays a stronger relationship with the baseline σ a ref (R 2 = .33) than the linear model (R 2 = .27). Figure 7 allows the identification of the best predictor for each local parameter. Given local parameters from a subset of plots, a linear relationship between them and their best predictor is derived and used to predict the value of the local In Subplots a and b, the black dots and dashed lines are used to illustrate the behavior of some plots as plotting all lines will make the graph unreadable. VCP, vertical coplanar mode parameters for the other plots. Those predicted local parameters are then used in one of the models (M1-M3). This process and the results are shown below for the multi-offsets (M2) and the multi-slopes (M3) models (M1 not shown). Hereafter, the subset of plots is composed of the 12 plots equipped with an ERT array as they are randomly distributed in the field. The choice of plots is somewhat arbitrary: another set of plots could have been selected, but they should span the largest possible range of σ and θ observed in the field.

Multi-offsets model
The multi-offsets (M2) model incorporates a local offset, b i , but a global slope, a g (Equation 9). As an illustration, Figure 8a compares, for a subset of plots (black line and dots), the multi-offsets model with its corresponding global model for VCP0.71. The global model compared here corresponds to Equation 6 where both slope, a g , and offset, b g , are uniform across the field. The multi-offsets model improves the accuracy of the predicted θ a compared with the global model (R 2 = .92 vs. .37) due to the inclusion of the local parameters b i (Figure 8a). Both models are fitted on all the plots available. In order to decrease the amount of data needed to obtain these local offsets, a linear relationship between the local offsets b i and the baseline θ a ref is derived using the data from a subset of plots (Figure 8b). This b i -θ a ref relationship is then used to predict b i for all the plots. Finally, in Figure 8c, those predicted offsets are used in the plocal multi-offsets model to obtain θ a . In this case, the R 2 of the multi-offsets model with the predicted parameters (.81) is better than for the global fit (.37).
The multi-offsets model focuses on the absolute values and not the differences. For the latter, the multi-slopes model is adapted further.

Multi-slopes model
The multi-slopes model (M3) presented in Figure 9 tries to fit a local model Δσ a and Δθ a (Equation 10). Figure 9a shows a comparison of the multi-slopes model and its global equivalent. In this case the global model contains a unique slope for the whole field. Similar to Figure 8, the introduction of a local parameter (slope a i ) improves the strength of the relationship from R 2 .71 to .86. In Figure 9b, a linear relationship is derived between the local slopes a i and the baseline σ a ref based on a subset of plots (R 2 = .64). This a i -σ a ref relationship is then used to predict the values of a i for all the other plots. Finally, those predicted slopes are used in Figure 9c in the multi-slopes model to predict Δθ a for all plots. The multislopes model with the predicted local parameters (plocal) has a higher R 2 (.68) than the global fit (.71). Figure 10 shows the quality of the prediction of M1, M2, and M3 using the predicted local parameters (plocal). The multi-offsets (M2) and multi-slopes (M3) models, which only have one local parameter, show better R 2 (M1: .16, M2: .53, and M3: .60) and a lower RMSE (M1: 0.04, M2: 0.02, and M3: 0.02) than the plocal linear model (M1), which has two local parameters. That means that the predicted soil moisture Vadose Zone Journal F I G U R E 9 Multi-slopes model fitted with differences in apparent values (VCP0.71). The gray dots show all the data available on the 216 plots.

Quality of the predicted local models
They represent the maximum number of information achievable if both electrical conductivity (σ) and soil moisture content (θ) are monitored on all the plots. In a more practical situation, only a subset of plots (black dots) are monitored for both σ and θ. Subplot a shows the multi-slopes model as well as a global relationship with a unique slope for all 216 plots (global). Subplot b shows the local slopes according to the baseline apparent soil electrical conductivity (σ a ref ). The black line corresponds to a linear relationship fitted on a subset of plots. This relationship is used to predict the local slopes for all the other plots. Subplot c shows the multi-slopes model using the predicted slopes from Subplot b (plocal). In Subplot a and b, the black dots and dashed lines are used to illustrate the behavior of some plots, as plotting all lines will make the graph unreadable. VCP, vertical coplanar mode F I G U R E 10 Quality of the predicted apparent soil moisture content (θ a ) vs. the observed θ a from (a) linear, (b) multi-offsets, and (c) multi-slopes models with predicted local parameters. The red line is the line of best fit with its 95% confidence interval (red shaded region). Both multi-offsets and multi-slopes models have one local parameter, whereas the linear model has two content from the multi-offsets (M2) or multi-slopes (M3) models is more accurate than from the linear model (M1).

Choice of the size of the subset of plots for plocal models
The size of the subset of plots needed for the plocal models needs to be chosen carefully. Figure 11 shows the effect of the number of selected plots on the RMSE of the prediction for the multi-offsets (M2) and the multi-slopes (M3) models. In this case, the RMSE does not change much if >10 plots are included in the subset.

Effect on the variety ranking
In a phenotyping context, we expect similarity in the rank of varieties whether observed (from neutron probe) or predicted (from EMI) soil moisture values are used. To assess the ranking improvement the predicted values of the global, local, and plocal models are averaged by variety. Then, the Spearman's rank correlation is computed between the observed and the predicted θ a (or Δθ a ). The Spearman's rank correlation has the advantage of being directly related to the ranking of the variety, which is a commonly used metric in crop breeding. A high value for this coefficient means, in our case, that higher predicted θ a is associated with higher observed θ a or that larger predicted θ a differences are associated with larger observed θ a differences, from examining absolute values or differences, respectively. Figure 12a shows the Spearman's rank correlations for the multi-offsets (M2) model using the baseline θ a ref as predictor of the local offsets. Figure 12b shows the Spearman's rank correlations for the multi-slopes (M3) model using the baseline σ a ref as predictor of the local slopes. Using the data in this study, the global models offer poor correlation compared to the local models, due to the Vadose Zone Journal F I G U R E 11 Effect of the size of the subset of plots on the predictions of the plocal (a) multi-offsets and (b) multi-slopes. After sorting the plots according to the baseline soil apparent electrical conductivity (σ a ), a subset of a given number of plots is selected at regular interval on the whole range of baseline values. θ obs and θ pred are observed and predicted soil moisture content. HCP, horizontal coplanar mode; VCP, vertical coplanar mode F I G U R E 12 Improvement in variety ranking in terms of the Spearman's rank correlation coefficient for (a) the multi-offsets and (b) the multi-slopes models. Each row of the table corresponds to a coil configuration. The columns are grouped by dates and subdivided into global, local and plocal models. The global models use field-specific parameters, and the local models use plot-specific parameters estimated using all the data available. The plocal model use the predicted plot-specific parameters estimated from baseline measurements (as in Figure 8b and Figure 9b). Bold numbers denote a significant correlation (p < .05). HCP, horizontal coplanar mode; VCP, vertical coplanar mode; θ a , apparent soil moisture content heterogeneity of the σ-θ relationship. This is true for all coil configurations. The plocal models (i.e., the models using the predicted local parameters) show higher correlation compared to their global equivalent. For the multi-offsets model, the improvement between global and the plocal is substantial (Figure 12a). When considering changes in soil moisture content (Figure 12b), the correlation with the global model is sometimes negative. This is a concern, as it means that an increase in σ a can be associated with a decrease in θ a after application of the global model. The local multislopes models increases this correlation substantially, especially for later dates. However, the plocal multi-slopes model shows relatively poor correlation even if it can compensate for the negative correlation observed in the global model in some cases.

Methodological limitations
The approach presented in this manuscript relies on apparent and not depth-specific electrical conductivity measurements to avoid the uncertainty arising from EMI inversion. Hence, we converted soil moisture content to apparent values using the practical cumulative sensitivity function (McNeill, 1980). However, the latter can have limitations especially on heterogeneous conductive soils. To estimate the errors that can arise from using the cumulative sensitivity function, Maxwell's equations can be used to reconstruct sensitivity functions based on a synthetic two-layer profile comparable with what is observed in the field (Callegary, Ferré, & Groom, 2007).
Both sensitivity functions are then used to compute the apparent soil moisture content. The maximum discrepancy between the two approaches is 0.01 cm 3 cm −3 , which is similar to the neutron probe accuracy (0.01 cm 3 cm −3 ). Given the magnitude of the errors, this probably has a more important impact on the changes in soil moisture content than on the absolute values. This might explain why the multi-slopes model works less well than the multi-offsets model in this study.
The dynamics of the soil moisture are complex, and isolating the effect of root activity is challenging. Whenever possible, measurements were collected at increasing potential soil moisture deficit and away from significant rainfall events (Figure 3). The drying observed in cropped plots compared with fallow plots suggests a substantial effect of the root activity ( Figure 3). However, the proposed approach does not aim at univocally measuring root water uptake, but rather at comparing soil moisture variation mainly induced by root activity for the different varieties.
The models described in the manuscript are simple linear models. More complex relationships can be used to relate soil moisture to electrical conductivity. For example, an exponential model was initially tested and showed similar performance to the linear model (see Section 3.5), hence the simplest model is chosen. In the linear models presented, the slope can be related to the soil surface conductivity, while the offset is more a function of the pore water conductivity. Both are functions of the soil texture and porosity (Garré et al., 2011). We do not have the information to investigate further the impact of these soil properties on the pedophysical parameters we derived for this field.
This study assumes that the samples taken on each plot (EMI, NP) are representative of the entire plot and that no substantial heterogeneity exists within the plot itself. Although we have no data to assess that this assumption is fulfilled for all the plots, the inverted ERT sections, which span 5.75 out of the 9 m of the plot length, suggest that this is the case.
One can question if the plot is the appropriate scale at which to investigate the variability of the θ-σ relationships. The use of variogram analysis can certainly help to determine the appropriate length scale at which the heterogeneity occurs. However, this method was not explored in this study, as our approach relies on the plot scale for practical reasons and to be consistent with additional phenotyping measurements at the site.
Finally, it has been assumed that the root system of the crop itself did not significantly contribute to the soil bulk apparent conductivity. Although there is evidence that suggests that coarser roots can affect the soil bulk electrical conductivity (Amato et al., 2008;Mary et al., 2017), finer herbaceous roots have been found to have a signal in magnitude similar to the effect of grain size or soil moisture content (Amato et al., 2009). Nevertheless, recent studies were able to isolate the electrical signature of roots themselves (Tsukanov & Schwartz, 2020). This could have great potential for phenotyping applications.

Ranking performance
Fitting a global model with field-specific parameters to all the data can lead to a satisfactory prediction of the soil moisture content, particularly if the differences expected between the treatments are large, such as for different types of vegetation (Jayawickreme et al., 2010), between fallow and cropped plots or between different soil types. However, when comparing a large number of similar varieties, this global model may be limited (Figure 12). In a phenotyping application, as here, using such a relationship may lead to false ranking of variates when using geophysical data ( Figure 12). As observed by Farahani, Buchleiter, and Brodahl (2005) for nonsaline soil, higher σ a is not always associated with greater soil moisture. Taking into account differences, it can also be seen that a large reduction in σ a is also not always associated with a large reduction in θ a . The negative correlations sometimes observed are of concern, as they lead to very different varieties ranking whether we consider σ a or θ a (Figure 12). The use of local parameters in the σ-θ relationship increases the Spearman's rank correlation for later dates as the soil moisture differences from the baseline become larger. The large number of parameters needed to fit the local models (linear, multioffsets, or multi-slopes) can be reasonably reduced using a relationship between the local parameters and the baseline σ a or θ a fitted on a subset of plots. The resulting plocal models that use those predicted parameters increase the accuracy of the prediction compared with global models (Figure 10). The R 2 value is often similar to or higher than those of the corresponding global models, but the ranking assessed (using the Spearman's rank correlation) is usually better (Figure 12). Note that the R 2 achieved are all below .6, which is relatively poor compared with what could potentially be achieved with a local relationship for all the plots (Figure 6b). Indeed, this improvement is mainly limited by the quality of the relationship between the local parameters and the predictors (Figures 8b and 9b). Hence, there is a need to select plots that span a wide range of conductivities to be monitored for both σ and θ (see Section 4.4) in order to have a more robust fit that is representative of the entire field.

Local models and parameters predictability
As seen in Figure 7, the offsets of the linear or multi-offsets models are mainly related to the baseline θ a . There is also a slight positive trend between the baseline σ a and the offsets of the linear model, but it is relatively weak compared with θ a , and it completely vanishes in the multi-offsets model. The The local slopes of the local linear model are well correlated with the baseline σ a . Considering differences, the multislopes model also shows good correlation between the local slopes and the baseline σ a . The conical shape of the data shown in Figure 4b and Figure 9a for the differences illustrates how different plots have different slopes. Stanley et al. (2014) show how the slopes of the σ-θ relationships vary between two sites with contrasting textures: sites with higher clay content, for example, result in greater values than those from sandier locations.
The multi-offsets and multi-slopes models have one contrasting assumption. The former assumes a unique slope for the entire field, whereas the latter uses plot-specific slopes. Having both plot-specific offsets and slopes leads to the local linear model, but its local parameters are difficult to predict using baseline measurement (Figure 7) and hence leads to poor estimates ( Figure 10). As the relationship between σ aθ a is largely offset dominated, we decided to fix the slope in the multi-offsets model to reduce the number of local parameters. For the differences, the effect of the offsets disappeared (Equation 3) and only the effect of the slopes has an impact on the relationship. This leads to the multi-slopes model (Equation 10).
As seen in Figure 13, the differences in observed σ a are well correlated with the baseline readings. Larger reductions in σ a are seen on plots with higher baseline σ a (Figure 13a).
Note that such a trend is not observed for θ a (Figure 13b). The fact that the σ a differences are still functions of the baseline reveals that the baseline σ a contains some information on how the σ a is likely to change: larger reductions are expected in areas of higher baseline σ a . This behavior explains why the starting σ a could be a good predictor of the slopes in the multislopes and linear models. Indeed, as the plots with higher baseline σ a show a larger increase in σ a with time for the same increase in θ a , they need to have a smaller slope to compensate. Smaller slopes are then found for higher baseline σ a (Figures 9b and 7). We believe this is related to the heterogeneity of the soil texture of the field where some areas are richer in clay than others.
Plots with higher baseline σ a tend also to have smaller offsets as well (Figure 7), but this relationship is not strong enough to be used for parameter prediction and θ a is preferred as the predictor (Figure 8b). Also, the prediction of the local parameters using the baseline readings is much better in the multi-offsets model (Figure 8b, R 2 = .82) than in the multi-slope model (Figure 9b R 2 = .64). This can explain why the multi-slopes model using predicted local parameters show only a slight improvement in variety ranking compared with the multi-offsets model ( Figure 12).
The multi-offsets and multi-slopes models are simplified ways to account for the variability due to the spatial heterogeneity of the θ-σ relationship. By reducing the number of local parameters compared with a local linear model, the local parameters are more correlated with baseline measurements and are thus easier to predict based on a subset of plots. In that way, they increase the ranking of the varieties and the accuracy of the predicted θ a compared with global models.

Improvement of the time-lapse approach
A key bottleneck in using the local models (M1-M3) is the predictability of the large number of local parameters they F I G U R E 14 Kernel density estimate (KDE) of the residuals for the (a) multi-offsets and (b) multi-slopes models for VCP0.71. For each, the global model represent a global (field-scale) linear relationship, whereas the local models use plot-specific parameters. The plocal model is the local model with the plot-specific parameters predicted from baseline apparent soil moisture content (θ a ) or soil apparent electrical conductivity (σ a ). VCP, vertical coplanar mode require. In this study, an approach was chosen where both variables (θ a and σ a ) are recorded on a subset of plots. In this case, the same 12 plots that served for the ERT calibration of the EMI data were arbitrarily chosen, as they are well distributed across the field and span the whole range of observed baseline values. In our case, a sample of 12 was large enough to reach the minimum RMSE achievable (Figure 11). Given the local parameters found on the selected plots, a relationship can be derived using the baseline σ a or θ a . This relationship can then be used to predict the values of the local parameters for the other plots. We believe that geostatistical tools can also be used to determine the number of sampling locations. However, we have not tested these in this paper.
Considering the above, we propose an improvement to the time-lapse approach described earlier to monitor the changes in soil moisture for large crop breeding experiment. After the first baseline EMI survey, plots with contrasting σ a are selected and equipped with soil moisture sensors (such as neutron probe access tube). The data collected on those plots will allow the estimation of the parameters for the multioffsets and multi-slopes models. Those parameters can then be expanded to the other plots using the baseline measurements (Figure 8b and Figure 9b).
The new approach is as follows: 1. Baseline survey on all the plots to acquire σ a ref and θ a ref : • multi-slopes: EMI with all configurations (σ a ref ) • multi-offsets: soil moisture measurements for all depths available to build an apparent soil moisture content measurements (θ a ref ) 2. Selection of plots with contrasting σ a to be equipped with θ sensors 3. Time-lapse EMI on all the plots and time-lapse θ on the selected plots: collection of multiple σ a -θ a datasets 4. Fit the multi-slopes (Equation 10) and multi-offsets (Equation 9) models on the selected plots to obtain the value of the local parameters: slope a i for multi-slopes and offset b i for multi-offset 5. Fit of linear relationship between those local parameters and the baseline value of the selected plots as in Figure 8b and Figure  This new approach offers a tradeoff between equipping all the plots with soil moisture sensors in order to fit a local models and using a unique global relationship for the entire field. Note that if a multi-offsets model is to be derived, baseline θ data are still needed, as they are the best predictors of the local offsets.

Analysis of the residuals
An increase in residuals can arise due the large number of local parameters. However, as Figure 14 shows, there is no substantial increase in the distribution of those residuals for the predicted local models compared with the global and local models. We can also see from Figure 8b and Figure 9b that even if the relationship is not perfectly fitted, the predicted parameters tend to stay in a reasonable range, avoiding the generation of outliers. Note that the residuals for the multislopes model are smaller than the residuals for the multioffsets, as the range of Δθ a (0 to −0.07) is smaller than the range of θ a (0.15 to 0.35).

CONCLUSIONS
High-throughput geophysical tools, in this case time-lapse EMI, offer great potential as a proxy measurement of soil moisture differences. When measurements are collected over increasing soil drying during crop growth, they may be linked to root activity in nonirrigated crop breeding field trials. The usual time-lapse approach is useful for removing the static effects of soil electrical conductivity but can be limited for ranking a large number of similar varieties in a heterogeneous environment. The spatial heterogeneity of the σ-θ relationship at the field scale has an impact on the ranking of the varieties and using a field-specific global relationship can lead to misleading interpretation. The proposed multi-offsets and multi-slopes models try to account for this heterogeneity by using plot-specific parameters that can be estimated from the baseline measurements. This improves the variety ranking between EMI and neutron probe data. A practical approach is proposed for such studies in which a baseline EMI survey is used to target sites for soil moisture monitoring, thus enhancing the ability to formulate predictions of the local σθ relationships. Although all the processing presented here was done with apparent conductivity measurements, the same process can be applied to depth-specific (inverted) measurements.