Simulation of spatial variability in crop leaf area index and yield using agroecosystem modeling and geophysicsbased quantitative soil information

Correspondence Cosimo Brogi, Agrosphere (IBG-3), Institute of Bioand Geosciences, Forschungszentrum Jülich, 52425 Jülich, Germany. Email: c.brogi@fz-juelich.de Abstract Agroecosystem models that simulate crop growth as a function of weather conditions and soil characteristics are among the most promising tools for improving crop yield and achieving more sustainable agricultural production systems. This study aims at using spatially distributed crop growth simulations to investigate how field-scale patterns in soil properties obtained using geophysical mapping affect the spatial variability of soil water content dynamics and growth of crops at the square kilometer scale. For this, a geophysics-based soil map was intersected with land use information. Soil hydraulic parameters were calculated using pedotransfer functions. Simulations of soil water content dynamics performed with the agroecosystem model AgroC were compared with soil water content measured at two locations, resulting in RMSE of 0.032 and of 0.056 cm3 cm−3, respectively. The AgroC model was then used to simulate the growth of sugar beet (Beta vulgaris L.), silage maize (Zea mays L.), potato (Solanum tuberosum L.), winter wheat (Triticum aestivum L.), winter barley (Hordeum vulgare L.), and winter rapeseed (Brassica napus L.) in the 1by 1-km study area. It was found that the simulated leaf area index (LAI) was affected by the magnitude of simulated water stress, which was a function of both the crop type and soil characteristics. Simulated LAI was generally consistent with the observed LAI calculated from normalized difference vegetation index (LAINDVI) obtained from RapidEye satellite data. Finally, maps of simulated agricultural yield were produced for four crops, and it was found that simulated yield matched well with actual harvest data and literature values. Therefore, it was concluded that the information obtained from geophysics-based soil mapping was valuable for practical agricultural applications.


INTRODUCTION
In recent years, significant improvements have been made in data acquisition, management, and modeling to support more productive and sustainable agriculture, which is vital to meet present food security challenges (Antle et al., 2017). Reductions in crop performance can be directly or indirectly attributed to a number of factors, such as insufficient nutrient availability, inadequate crop and soil management practices, adverse weather conditions, occurrence of pests and diseases, and water shortage (Baret, Houles, & Guérif, 2007;Foley et al., 2011;Sánchez, 2010;Slingo, Challinor, Hoskins, & Wheeler, 2005). Worldwide, agriculture is the largest consumer of freshwater (Brauman, Siebert, & Foley, 2013), with most crops being produced in rainfed conditions (Rosegrant, Ringler, & Zhu, 2009). Within this context, water stress caused by below-average precipitation is considered to be a severe threat for sustainable crop production in the coming decades (Anjum et al., 2011). It will therefore be essential to produce more food per unit of water (Brauman et al., 2013). Unfortunately, the occurrence of water stress in crops cannot always be prevented, especially if irrigation infrastructure is not available (Ceglar, Toreti, Lecerf, Van der Velde, & Dentener, 2016;Prasad, Staggenborg, & Ristic, 2008). For most crops, water stress results in a decrease of leaf production, a reduction in leaf area index (LAI), an increase in senescence rate, and finally yield reduction (Baret et al., 2007). In the case of water stress, actual crop growth and yield are influenced by factors such as the amount and timing of precipitation, soil physical properties, and soil horizonation, since all these factors have an effect on the plant-available soil water content (Krüger et al., 2013;Lück, Gebbers, Ruehlmann, & Spangenberg, 2009). Nevertheless, our ability to simulate within-field variability of water stress and the resulting reduction in crop growth is still limited (Batchelor, Basso, & Paz, 2002).
One effective method for estimating crop water stress is the use of process-oriented crop growth models (Batchelor et al., 2002). In this approach, the interaction of water stress with plant growth is simulated with a high temporal resolution, thus enabling yield predictions in various environments and meteorological conditions. Among others, such models have been successfully applied to determine the effects of the spatial variability of soil moisture on yield Paz et al., 1998Paz et al., , 1999, to investigate yield loss in corn (Zea mays L.) in a cool climate due to abiotic stress (Žydelis, Weihermüller, Herbst, Klosterhalfen, & Lazauskas, 2018), and to validate management zones derived from satellite-based vegetation indices (Basso, Ritchie, Pierce, Braga, & Jones, 2001). These process-oriented crop growth models typically rely on a one-dimensional description of water flow in the soil column  and require a detailed description of the soil profile characteristics including soil hydraulic properties (Boenecke, Lueck, Ruehlmann, Gruendling, & Franko, 2018). In general, this information is obtained from generalpurpose maps (Boenecke et al., 2018) that discretize soil in relatively large polygons. Typically, such maps are considered to be too coarse for precision agriculture and within-field Core Ideas • Quantitative information on soil improves largescale agricultural simulations. • Water stress caused by soil characteristics affected crop simulations. • Water stress decreased growth and yield of the simulated crops. • Detailed thematic agricultural maps are relevant in practical agricultural applications.
management (Robert, 1993) because site-specific soil heterogeneity is not adequately captured (Gebbers & Adamchuk, 2010;Heuvelink & Webster, 2001). As a result, a more accurate spatial description of soil properties is frequently needed (Krüger et al., 2013) to perform effective simulations within an agricultural context. A possible strategy to address the need for high-resolution soil information is geophysics-based soil mapping. It has been shown that a range of geophysical properties provide useful proxies for soil properties in precision agriculture (Adamchuk, Hummel, Morgan, & Upadhyaya, 2004;Allred, Daniels, & Ehsani, 2008;Gebbers & Lück, 2005;Grisso, Alley, Holshouser, & Thomason, 2005;Vitharana, Van Meirvenne, Simpson, Cockx, & De Baerdemaeker, 2008). Geophysical measurements have been particularly useful in the definition of management zones (Brogi et al., 2019;Galambošová, Rataj, Prokeinová, & Presinska, 2014;King et al., 2005;Oldoni & Bassoi, 2016;Taylor, Wood, Earl, & Godwin, 2003). Krüger et al. (2013) used electromagnetic induction (EMI) and ground-penetrating radar (GPR) to characterize site-specific variations in soil properties, which improved the simulations of soil water dynamics and biomass production on a 4.4-ha field. Wong and Asseng (2006) used EMI to map the plant-available water content and used simulations to illustrate how variations in available water interacted with the amount and timing of precipitation and yield variability within a single 70-ha field. Similarly, Boenecke et al. (2018) used EMI as a basis for simulating the spatial variability of soil water content and yield at the farm scale (30 ha). Despite these successful examples, there is a general lack of studies that link soil maps with modeling applications at scales larger than a single farm (Krüger et al., 2013) and for multiple crops. Successful simulations of the spatial variability of water stress and the associated decrease in crop productivity within large and complex agricultural environments could lead to significant improvements in our knowledge of the soil-plant system, thus providing valuable support in developing long-term strategies for decisions in agricultural management (Krüger et al., 2013).

Vadose Zone Journal
Within this context, the aims of this study are (a) to set up and validate a spatially distributed crop growth model for an agricultural area of 1 × 1 km 2 composed of 51 fields cultivated with different crops using the high-resolution geophysicsbased soil map obtained by Brogi et al. (2019), and (b) to spatially assess simulated yield relative to the yield obtained with optimal water supply for the given climatic conditions. For this purpose, the geophysics-based soil map was intersected with land use information to obtain 80 unique soilcrop combinations. Soil hydraulic parameters were estimated using pedotransfer functions (PTFs) and used to run the agroecosystem model AgroC to simulate vertical water fluxes and crop growth for each unique soil-crop combination. Simulations of soil water content were compared with measurements obtained at two locations within the study area to assess the reliability of the model parameterization. In a next step, AgroC simulations of the temporal development of the LAI of sugar beet (Beta vulgaris L.), silage maize, potato (Solanum tuberosum L.), winter wheat (Triticum aestivum L.), winter barley (Hordeum vulgare L.), and winter rapeseed (Brassica napus L.) were compared with LAI calculated from normalized difference vegetation index LAI NDVI that was determined from RapidEye multispectral satellite data acquired between April and September 2016. Finally, maps of simulated agricultural yield were produced to support decision making in agricultural management.

Study area
The study area is located in the Rur catchment near Selhausen, Germany (50 • 51 ′ 56 ″ N 6 • 27 ′ 03 ″ E). It is a square area of ∼1 × 1 km ( Figure 1). The average annual precipitation is 715 mm and the mean annual temperature is 10.2 • C . The area is composed of 50 fields cultivated in rotation with silage maize, potato, winter rape, winter wheat, winter barley, and sugar beet. Occasionally, oat (Avena sativa L.) is cultivated and some fields are left bare or are covered with grass. The area is characterized by a heterogeneous agricultural management, since >20 different landowners are managing different crop rotations (Brogi et al., 2019).
The shallow geology of the area consists of two main features: the eastern upper terrace and the western lower terrace of a paleo-river system. The upper terrace has an altitude of ∼110-113 m asl, whereas the lower terrace has an altitude of ∼101-103 m asl. The two terraces are separated by a 10% slope that dips westbound and is oriented north-northwest to south-southeast. Both terraces are characterized by quaternary sediments. The upper terrace is composed of Pleistocene sand and gravels (Röhrig, 1996) buried by aeolian sediments of various thickness (Klostermann, 1992;Patzold et al., 2008; F I G U R E 1 Satellite image of the study area (DigitalGlobe: September 2016) showing the Universal Transverse Mercator coordinates (UTM), the fields that are included in the study and the locations of the in situ leaf area index (LAI) measurements (green dots), the water table measurements, the meteorological station, and the soil water content measurements (points P01 and P02) Vandenberghe & Van Overmeeren, 1999). The lower terrace is composed of Pleistocene loess and translocated loess from the Holocene. Locally, sand and gravels from the Pleistocene and Holocene are found within 2-m depth below the loess sediments (Brogi et al., 2019;Röhrig, 1996). The dominant reference soil groups are Cambisols, Luvisols, Planosols, and Stagnosols (WRB, 2015). Previous research in this area showed that crop performance during long periods of water scarcity is strongly influenced by soil heterogeneity at the field scale and beyond (Brogi et al., 2019;Rudolph et al., 2015;Stadler et al., 2015;von Hebel et al., 2018). The study area is part of the Terrestrial Environmental Observatories (TERENO) network (Bogena et al., 2018;Schmidt, Reichenau, Fiener, & Schneider, 2012;Simmer et al., 2015). Table 1 provides information on the available meteorological measurements performed in field F11 (measurement location shown in Figure 1).
Within the study area, a 2.32-ha field cropped with sugar beet and a 9.54-ha field cropped with barley were investigated in more detail (fields F01 and F11 in Figure 1). Crop yield of both fields in 2016 was provided by the respective farmer. In field F01, the available fresh weight of harvested  Figure 1) to monitor volumetric soil water content from 28 Apr. 2016 to the 18 Oct. 2016. At each location, two SMT-100 soil water content sensors (Truebner) were installed at three depths (10, 20, and 50 cm). These sensors were calibrated to provide soil dielectric permittivity before installation (Bogena, Huisman, Schilling, Weuthen, & Vereecken, 2017), and the volumetric soil water content was obtained using the equation of Topp, Davis, and Annan (1980).

Geophysics-based soil map
A geophysics-based soil map that provides detailed soil information of the study area ( Figure 2a) was generated by Brogi et al. (2019). This high-resolution soil map was obtained using multiconfiguration EMI measurements made with a CMD Mini Explorer, Special Edition (GF Instruments) directly after harvest of each field in 2016. The EMI survey resulted in six maps of the apparent electrical conductivity (ECa), each with a different depth of investigation. First, the six ECa maps with a resolution of 1 m 2 were stacked in a multiband image. In a following step, 18 soil units were identified using information contained in the ECa maps, commonly available soil maps, and expert knowledge from previous studies and field observations. Following a field-by-field workflow, areas belonging to a specific soil class were identified in each field (so-called training areas) and a maximum likelihood was used to assign each pixel to the most statistically probable class. The application of this supervised classification methodology for each field resulted in a high-resolution map with 18 areas with similar EMI response. In a next step, 100 sampling locations were visited in January 2017. At each location, horizon type and thickness were recorded. In addition, a total 552 soil samples were collected for different locations and horizons to obtain information on the grain size distribution. Finally, all profile descriptions and grain size distributions were averaged to obtain a representative soil profile with quantitative information on soil texture for each of the 18 soil units. A good match was found when these soil units were compared with patterns in crop stress visible from satellite in 2015. The 18 soil units are divided in four groups. Group A (soil units A1a-d) represents units located in the upper terrace where aeolian sediments of various thickness are found on top of coarse sand and gravel material. Group B (soul units B1a-b and B2a-c) represents units located on the slope that divides the upper and the lower terrace. Group C (soil units C1a-b and C2a-b) is composed of soils of the lower terrace where only aeolian sediments are generally found in the top 2 m. These sediments are locally buried with anthropogenic material (units C2a-b). Group D (soil units D1a-d and D2a) represents units of the lower terrace where aeolian sediments of various thickness are found on top of coarse sand and gravel material. Locally, anthropogenic material is found within the soil profile (unit D2a). For more information on the derivation and validation of this high-resolution soil map based on geophysical measurements, the reader is referred to Brogi et al. (2019).

Land use and unique soil-crop combinations
Information on crop type, emergence, and harvest dates was recorded during 2016. The land use map shown in Figure 2b was produced by combining field geometries retrieved from satellite images (ESRI, 2018) with in situ differential global positioning system (DGPS) measurements (Trimble). The crops present in 2016 were sugar beet, silage maize, potato, winter wheat, winter barley, and winter rapeseed. The total area for each crop is provided in Table 2. A small percentage of the area was covered with bare soil, oat, and grass. These fields were excluded from analysis. Different emergence and harvest dates were recorded in separate fields for the same  crops. However, these differences were rather small (e.g., 3-10 d for emergence and 1-5 d for harvest dates). Therefore, it was assumed that each crop type is characterized by a single harvest date and a maximum of two different emergence dates in our simulations (see Table 2). In a next step, the soil and land use map were intersected. This resulted in 80 unique soil-crop combinations (Figure 2c).

Estimation of hydraulic parameters
The PTFs provided by Rawls and Brakensiek (1985) were used to estimate the hydraulic parameters of each layer of the soil profiles provided by the geophysics-based soil map. This PTF uses the fractions of sand, silt, and clay, and the dry bulk density to estimate the hydraulic parameters. The dry bulk density of the fine fraction <2 mm (BD <2 ) was not directly determined during soil profile characterization and sampling for soil texture. In the upper plow layers, a lower BD <2 of 1.30 g cm −3 for Ap horizons and a BD <2 of 1.40 g cm −3 for AB horizons were assumed due to regular tillage (Ehlers, Köpke, Hesse, & Böhm, 1983;Unger & Jones, 1998). For the deeper soil horizons, a BD <2 of 1.50 g cm −3 was assumed for fine sediments, and a BD <2 of 1.60 g cm −3 was assumed for 2C horizons that contained a lot of coarse material. These assumptions were based on results from previous sampling campaigns conducted on two fields within the study area (fields F01 and F11 in Figure 1). Because of the high gravel content of some soils, the bulk density of the fine earth fraction was corrected for gravel content according to Brakensiek and Rawls (1994) by where BD t is the bulk density of the soil, BD >2 is the bulk density of gravel, and G v is the volume of gravel (Flint & Childs, 1984). The bulk density of gravel was assumed to be that of solid quartz: BD >2 = 2.65 g cm −3 (Brakensiek & Rawls, 1994). The volume of gravel was calculated using where G w is the weight fraction of gravel (Flint & Childs, 1984), which was available from the soil texture information provided by the geophysics-based soil map. Due to the high gravel content of the 2C horizon (estimated to be ∼70%), the estimates of the hydraulic properties were deemed to be less accurate for this horizon. To avoid introducing variation into the simulation results due to these uncertain hydraulic properties, it was decided to set the hydraulic parameters of the 2C horizon to the same values for all soils of the upper terrace (soil units A1a-d) and for all soils of the lower terrace (soil units D1a-d and D2a). In both cases, the hydraulic parameters were estimated from the average soil texture and corrected bulk density obtained from the estimation procedure outlined above. Finally, the estimated saturated hydraulic conductivity (K s , cm h −1 ) of the 2C horizon was also corrected for gravel content according to Brakensiek and Rawls (1994) were K b is the saturated hydraulic conductivity of the bulk soil (fine earth and gravel) and K s is the saturated conductivity of the fine earth fraction estimated by the PTF of Rawls and Brakensiek (1985). An example of the estimated hydraulic parameters is shown in Table 3 for soil unit A1a. For comparison, the hydraulic parameters of soil unit A1a obtained by applying the same procedure but using the commonly applied ROSETTA PTF (Zhang & Schaap, 2017) are also shown and are used later on for comparison.

Estimation of leaf area index from satellite observations
Leaf area index was estimated from RapidEye multispectral satellite imagery using a logarithmic relationship to calculate LAI NDVI from vegetation indices . Successful applications of this procedure have been published by Hasan et al. (2014), Montzka et al. (2016), Post et al. (2018), Reichenau et al. (2016), and Rudolph et al. (2015). For this study, six RapidEye Level 3A images provided with radiometric, sensor, and geometric correction and with a resolution of 5 m were available. The six images were acquired on the 14 Mar., 20 Apr., 28 May, 9 June, 12 Aug., and 8 Sept. 2016 and thus cover the full growing season of various crops. Further data had to be discarded because of dense cloud cover and relatively poor illumination conditions during winter with associated low signal-to-noise ratios.
In a first step, the normalized difference vegetation index (NDVI) was calculated for each image pixel. Then, the fractional vegetation cover (FVC NDVI ) was calculated from NDVI using the method proposed by Ali et al. (2015). A detailed description of this process is presented in Appendix A. Finally, the LAI NDVI was calculated using where k(θ) is the light extinction coefficient for a given solar zenith angle (θ). This coefficient is a measure of attenuation of radiation in the canopy and depends on factors such as the latitude, date, solar elevation and declination, terrain geometry, as well as canopy structure and architecture (Campbell, 1986;Norman & Campbell, 1989;Propastin & Erasmi, 2010;Ross, 2012). The k(θ) is challenging to measure in the field and is generally assumed to be constant in most biogeochemical models due to a lack of in situ measurements (Zhang, Hu, Fan, Zhou, & Tang, 2014). For this reason, k(θ) is generally obtained from literature values or retrieved inversely by comparing satellite observations with in situ LAI measurements (Propastin & Erasmi, 2010). In this study, k(θ) was calibrated to match RapidEye observations with in situ destructive LAI measurements of winter wheat, winter barley, and sugar beet collected throughout the growing season. For this, a set of 45 in situ destructive LAI measurements obtained between 22 Mar. and 7 Sept. 2016 was used. The fields in which these measurements were performed are shown in Figure 1a. At each sampling location, destructive measurements of LAI were performed by sampling 0.5 m of plants from three adjacent rows in wheat and barley. For sugar beet, three individual plants were collected and crop density was determined in the field. Subsequently, the leaf area was measured in the laboratory by using a flatbed scanner (Epson GT-15000, Seiko Epson Corporation) and a publicdomain image analysis software (ImageJ, National Institute of Health). Upscaling to a square meter was performed based on the harvested area. These in situ LAI values were compared with the respective FVC NDVI values to obtain an estimate of k(θ) used in the calculation of LAI NDVI . In this process, only in situ measurements that were performed 8 d or less before or after the acquisition date of a satellite record were used. The number of in situ measurements, the R 2 between FVC NDVI and in situ LAI, and the k(θ) of each crop are shown in Table 4.
T A B L E 3 Hydraulic parameters of soil unit A1d estimated using the Rawls and Brakensiek (Ra. & Br.) pedotransfer function (PTF) (Rawls & Brakensiek, 1985a)  In sugar beet, the initial R 2 between in situ measurements and calculated LAI NDVI after calibrating k(θ) was rather low at .29. Therefore, in this study, two different values of k(θ) were used to optimize the correlation between LAI NDVI and in situ LAI: one value for RapidEye data from March to August, and one value for September (Table 4). This improved the resulting R 2 considerably to .79 (May-August) and .74 (September). The use of two different k(θ) values within the growing season can be justified by (a) the fact that the relationship between NDVI and LAI saturates at high LAI values once leaves completely cover the soil towards the end of the growing season, and (b) a change in the geometry and in the greenness of the leaves, and therefore in NDVI, that can occur towards the end of the growing season and/or during water scarcity. In fact, nonirrigated sugar beet generally suffers from water stress (Hoffmann & Kenter, 2018), which is known to have an immediate effect on leaf greenness, geometry, and other plant characteristics (De Costa & Dennett, 1992) affecting NDVI and therefore LAI NDVI . Late-summer water stress in sugar beet was reported in previous studies within the study area (Brogi et al., 2019;Rudolph et al., 2015;Stadler et al., 2015). Similarly, two k(θ) values were used in winter wheat, since the R 2 obtained using all available measurements was only .6. When two different values of k(θ) were used in March-April and in May-July, much improved R 2 values of .90 and .98 were obtained. On the contrary, barley showed an overall higher R 2 of .75 when all available LAI measurements were used to calibrate k(θ). In the case of maize, potato, and winter rapeseed, k(θ) was set to .25, since no ground-based LAI measurements were available for calibration. This value was suggested by Ali et al. (2015) to obtain LAI NDVI maps from RapidEye satellite data in this area.

The AgroC model and simulation setup
AgroC is an agroecosystem model that couples (a) the SOILCO2 module (Šimůnek & Suarez, 1993;Šimůnek, Suarez, & Šejna, 1996) for simulating vertical water, heat, and CO 2 fluxes in a soil column, (b) the RothC module for simulating the turnover of organic C (Coleman & Jenkinson, 1996), and (c) the SUCROS module for simulating crop growth and organ-specific dry matter (DM) accumulation in crops (Spitters, van Keulen, & van Kraalingen, 1989). The coupling and validation of the SOILCO2 and the RothC module is reported by Herbst et al. (2008). The development and validation of the AgroC model is described in detail in Klosterhalfen et al. (2017). A short overview of the key components of the AgroC model relevant for this work is provided in Appendix B.
In the SUCROS module, a root water uptake stress factor φ(h) reduces the potential root water uptake in dry and wet soil conditions. This reduction factor is calculated according to Feddes, Kowalik, and Zaradny (1978): where h 0 , h 1 , h 2 , and h 3 (cm) are prescribed threshold pressure heads. In this study, the threshold pressure head values were set to h 0 = 0, h 1 = −20, h 2 = −5,000, and h 3 = −16,000 cm according to Vanclooster, Viaene, and Diels (1995). The simulation of root water uptake stress affects estimated C assimilation and biomass production as described in detail in Appendix B. The soil column for the AgroC simulations was discretized with a maximum node spacing of 1 cm for the deeper soil and a finer spacing near the soil surface (0.1 cm). In general, two different setups were used. The first setup was used for soil profiles that have a coarse 2C horizon at the bottom of the soil profile. This includes all soils of the upper terrace and some soils of the lower terrace (i.e., soil units A1a-d, D1a-d, and D2a in Figure 1). This coarse and less conductive layer is known to have a strong influence on water flow and crop growth (Brogi et al., 2019;von Hebel et al., 2018). For these soils, the depth of the simulation domain was defined to be from the surface down to a depth of 3 cm into the C2 horizon. This resulted in different depths for the simulation domains that ranged from 52 to 137 cm depending on the soil profile description. The lower boundary condition was set to free drainage for these profiles. This setup and selection of boundary conditions was required to avoid unrealistic upward water flow from the C2 horizon that occurred for deeper soil profiles. The second setup was used for soil profiles where fine sediments were present at depth. In this case, the simulation domain extended from the surface down to a depth of 2.0 m. A water table with variable depth was used to define a variable pressure head as the lower boundary. The maximum and minimum water table depth was obtained from a CTD-5 sensor (Decagon Devices) installed in a well in the western part of field F10 (Figure 1). In 2015 and 2016, the water table depth varied periodically over the year. The minimum depth was 2.0 m on ∼15 January and the maximum depth was 2.6 m on ∼15 July.
Crop-specific input parameters required for AgroC were mainly obtained from literature values (Allen, Pereira, Raes, & Smith, 1998;Bolinder, Angers, & Dubuc, 1997;Boons-Prins, De Koning, & Van Diepen, 1993;Borg & Grimes, 1986; Penning de Vries, Jansen, ten Berge, & Bakema, 1989;Spitters et al., 1989;Van Heemst, 1988;Vanclooster et al., 1995). For each crop, a specific maximum rooting depth was set: 150 cm for sugar beet and silage maize, 140 cm for rapeseed, 120 cm for barley, and 100 cm for wheat and potato. The root distribution with depth was calculated using the method of Schamschula, Barmes, Keyes, & Gulbinat (1974). Due to the high bulk density of the 2C horizons, it is assumed that roots cannot grow into this layer due to high penetration resistance (Daddow & Warrington, 1983). Therefore, the rooting depth was reduced when a coarse 2C horizon was present within the rooting depth. In such a case, the root distribution between the surface and the 2C horizon was extracted from the distribution that was calculated with the crop-specific maximum rooting depth (i.e., a truncated root length distribution).
Meteorological data for 2015 and 2016 were used to define the upper atmospheric boundary condition using precipitation and calculated potential grass reference evapotranspiration according to the Penman-Monteith approach (Allen et al., 1998). All simulations started on 1 July 2015 and ended on the 31 Dec. 2016, thus allowing simulation of both summer and winter crops. Two strategies were again used to define the initial pressure heads within each soil profile. The first strategy was applied to soil profiles where a water table with variable depth was used as lower boundary condition. Here, a spin-up simulation was used by repeatedly running a period of two years (1 Jan. 2015 to 31 Dec. 2016) until the pressure head within the soil column did not change between subsequent simulations. This spin-up strategy resulted in very low water content within the soil column when simulating soil profiles with a free drainage boundary at the bottom. To avoid such unrealistically low initial water contents in these soil profiles, the pressure head of the underlying 2C horizon was set to −1 cm and hydrostatic equilibrium was assumed for the rest of the profile for initialization. Figure 3 shows the comparison between measured and simulated soil water content dynamics at locations P01 and P02 (see Figure 1). For this, horizonation and soil texture were obtained from the geophysics-based soil map, whereas bulk density was obtained from literature values and from previous unpublished soil sampling in the area. All these data were than propagated through two PTFs to obtain the soil hydraulic parameters used in the AgroC simulation. In Figures 3a-3b, the soil water content measured at P01 and P02 is compared with the soil water content for soil units A1a and A1d simulated with AgroC by using the PTF provided by Rawls and Brakensiek (1985). In general, the simulated and measured soil water content showed a similar response to atmospheric forcing in both locations with increasing water content after precipitation events followed by a dry-down in periods without precipitation. For the location P01 in soil type A1a (Figure 3a), the average RMSE between measured Vadose Zone Journal F I G U R E 3 (a) Simulated and measured soil water contents (WC, cm 3 cm −3 ) for sugar beet grown in soil unit A1a at P01, (b) simulated and measured soil water contents (cm 3 cm −3 ) for sugar beet grown in soil unit A1d at P02 based on the hydraulic parameters estimated by the pedotransfer function (PTF) of Rawls and Brakensiek (1985b), simulation performed for sugar beet grown (c) in unit A1a at P01, and (d) in unit A1d and P02 based on the ROSETTA PTF and simulated water content for all three soil depths was 0.056 cm 3 cm −3 . In general, the simulated water content was very similar for all three depths, whereas the measured water content considerably increased with depth and showed a stronger response to atmospheric forcing. For location P02 in soil type A1d (Figure 3b), the average RMSE between measured and simulated water content for all three soil depths was considerably lower at 0.032 cm 3 cm −3 . However, differences between measured and simulated water content were also clearly present. For example, two peaks in measured soil water content at 10-cm depth in August were not well captured by the model. Also, the measured soil water content at 50 cm was overestimated during the dry period from August onward.

Water content simulations for field F01
This difference in performance for the two locations is related to discrepancies between the actual soil profile and the average soil profile of the two soil units obtained from the geophysics-based soil map. In general, the actual horizon depths at point P02 were similar to the horizon depths of the soil profile of soil unit A1d. At this location, a rather thin layer of loess sediment was deposited on a coarse sand and gravel layer. The actual depth of the loess layer was 47 cm, whereas it was 49 cm in soil unit A1d of the soil map. Therefore, there was a good match between measured and simulated water content. For point P01 in soil unit A1a, a larger mismatch between the actual and the average soil profile was observed. In particular, the actual thickness the loess layer at P01 was 160 cm, whereas the average thickness of this layer was only 87 cm in soil unit A1a. The differences in soil profile description are a likely explanation for the mismatch between the simulated and the measured water content at this location. Obviously, these differences are due to the scale mismatch between the point measurements and the geophysics-based soil map. This soil map is expected to capture the main variability in soil properties at the 1-km 2 scale but is not expected to capture smallscale variabilities in profile depth and description within one soil unit. Apparently, P01 was not a representative location for the soil unit A1a, and therefore it could not be simulated with the same accuracy as P02.
To investigate the sensitivity of the AgroC simulations to the choice of the PTF, the simulations for locations P01 and P02 were repeated with soil hydraulic parameters estimated from the widely used ROSETTA PTF (Figures 3c-3d). This resulted in a stronger mismatch between measured and simulated water content. In particular, the RMSE increased from 0.056 to 0.134 cm 3 cm −3 for soil type A1a, and from 0.032 to 0.072 cm 3 cm −3 for soil type A1d. Also, the simulated water content obtained with the ROSETTA PTF seemed unrealistically low for both soil units. The estimated hydraulic parameters of soil unit A1d using the Rawls and Brakensiek as well as the ROSETTA PTF are given in Table 3. As expected, all soil hydraulic parameters showed differences when estimated by the two PTFs (Table 3), whereby the largest differences were obtained for K s (up to five times higher K s estimated by the ROSETTA PTF). This higher K s likely explains the low simulated water contents based on the ROSETTA PTF and therefore the increase in the mismatch between measurements and simulation. Since the unrealistically low simulated water content was observed for all soil units, it was preferred to use the Rawls and Brakensiek PTF in this study. Overall, it is concluded that the presented parameterization strategy provides reasonable predictions of soil water content dynamics considering that the soil profiles were obtained from a soil map and no calibration was used to improve the fit between measured and modeled soil water content. Alternative methodologies could be implemented in future research to improve automatization in the estimation the soil hydraulic parameters (e.g., by using inversion strategies). However, such a methodology was not pursued in this large-scale study since continuous measurements of soil water content were only available at two locations and were representative for only two soil units, thus preventing the parametrization of the other 16 units of the geophysics-based soil map using model inversion.

Validation of leaf area index and yield simulations
In a next step, the results of the AgroC simulations at the square kilometre scale were evaluated using LAI NDVI derived from satellite remote sensing. When available, yield information was also used for model evaluation. The results will be presented for three groups of units: (a) the soils of the upper terrace (units A1a-d), (b) the soils of the lower terrace with underlying fine layers (soil units B1a-b, B2a-c, C1ab, and C2a-b), and (c) the soils of the lower terrace with underlying coarse layers (soil units D1a-d and D2a). For simplicity, we will refer to the second group as BC soil units. Below, we present results for the summer crops that are generally more prone to water stress (sugar beet, silage maize, and potato), followed by the winter crops (wheat, barley, and rapeseeds).

Simulation of sugar beet
Figures 4a-4c shows the simulated and observed LAI NDVI for all soil units grown with sugar beet. The observed LAI NDVI was obtained by averaging all LAI NDVI values within each unique soil-crop combination. High LAI NDVI was observed for all BC soil units, which indicates that the crops in this area did not suffer from water stress. A reduction in LAI NDVI was observed in soil units A1a-d and D1a-d, and the reduction appears to be proportional to the depth of the coarse layer. In order to capture these differences in LAI between the three groups of soil units with the AgroC model, it was necessary to use three different crop parameterizations. In fact, the water stress factor simulated with AgroC reduced the accumulation of new dry mass but had less influence on the simulated green LAI, which was accumulated before the occurrence of prolonged periods of water stress in July and August (Figure 4d). Additional reduction in LAI NDVI during periods of water stress might be caused by curling and wilting of leaves, an aspect that is also not implemented in AgroC. In order to match the differences in observed LAI between the three groups of soil units, the death rate of the leaves was calibrated by increasing this rate in soil units A1a-d and D1a-d until a meaningful drop in simulated LAI in August-September was obtained. In the final set of simulations, three different rates were used (one in each group of soil units). The death rate of the leaves was highest in soil units A1a-d, intermediate in D1a-d, and lowest in the BC units (see Appendix C for a detailed description of the input parameters). Since the death rate of the leaves was the same within each group of soil units and all other plant parameters were identical in all simulations, the differences in LAI shown within each group of soil units are solely due to differences in soil parameterization in terms of layering and hydraulic properties.
In soil units A1a-d (Figure 4a), the LAI NDVI was similar for all four soil units early in the growing season (May-June). This was well reproduced by the model simulations, despite the strong water stress factor in this period (Figure 4d). Apparently, this stress period did not strongly affect the simulated accumulation of aboveground biomass at this initial stage of crop development. In the second half of June and July, there was a period with no water stress where the simulated LAI increased rapidly. Afterwards, the absence of rain resulted in water stress and an associated reduction of LAI in August and September. Here, the simulated LAI was lower in soil units with higher water stress and matched well with both the maximum value and the decrease of observed LAI NDVI in response to water stress.
In soil units BC (Figure 4b), simulated LAI and observed LAI NDVI generally matched well, but only showed limited variability between simulations of different soil units. This is explained by the lack of simulated water stress in this Vadose Zone Journal

F I G U R E 4 Simulations for sugar beet. Simulated leaf area index (LAI, lines) and observed LAI calculated from normalized difference vegetation index (LAI NDVI , points) (a) in soil units A1a-d, (b) in soil units BC, and (c) in soil units D1a-d (confidence interval of LAI NDVI shown with whiskers bars). Simulated water stress factor (d) in soil units A1a-d, and (f) in soil units D1a-d.
Note that a water stress factor of 1 does not produce a reduction of biomass accumulation (no water stress), whereas a factor of 0 represents the strongest simulated water stress. (e) Simulated mass of storage organs in soil units A1a-d and BC. Letters along the x axis represent months of the year, from March to November. Non-str., nonstressed area. A mismatch can be observed on 9 June where simulations matched the LAI NDVI of soil units C1a-b and C2a-b (LAI = 1.1-1.4) but did not match the LAI NDVI of the soil units B2a-b and B2a-c (LAI = 2.2-2.6). This is attributed to differences in seeding and emergence dates between soil units that are not captured in the model. This is supported by the fact that the two groups of soil units are located in separate fields with different owners (units B1a-b or B2a-c in field F47 and units C1a-b or C2a-b in fields F12 and F50).
In soil units D1a-d (Figure 4c), a general decrease in LAI proportional to the intensity of simulated water stress factor was observed in the second half of the growing season (similar to the soil units A1a-d). However, the simulated stress in D1a-d is generally lower than in A1a-d since the coarse layers are generally deeper in D1a-d and have different soil hydraulic parameters (Figure 4f). Again, simulated LAI and observed LAI NDVI matched well. Until August, the simulated LAI for soil units D1a-c was very similar. A considerable difference between soil units was simulated on 8 September because water stress reduced the LAI of soil units D1a-c. This simulated behavior is well supported by the observed LAI NDVI . On the contrary, the simulations for soil unit D1d underestimated LAI NDVI , likely because water stress was overestimated. However, this underestimation occurred in a small portion of the total area cropped with sugar beet in 2016, since soil unit D1d represented only 0.3 ha from a total of 26.5 ha.
The simulated dry weight of storage organs of sugar beet expressed in tons per hectare is shown in Figure 4e for selected soil units. The simulated weight of storage organs at harvest of the soil units A1a-d ranged from 12.0 to 15.9 t ha −1 and was well below the simulated weight of storage organs of soil units BC (20.7 t ha −1 ). Here, a reduction in the mass allocated to the storage organs is apparent from July until the end of the growing season, which was proportional to the magnitude of water stress. The simulated yield of soil units A1a-d was compared with the actual yield at harvest of field F01 in 2016. The area-weighted average simulated yield for soil units A1ad in field F01 was 14.3 t ha −1 and well matched the actual yield of 14.2 t ha −1 . The simulated weight of storage organs at harvest of soil units BC and D1a-d was compared with the average yield at harvest of field F11 in 2011, 2014, and 2017. Again, the average simulated yield of 19.7 t ha −1 in 2016 was similar to the average observed yield of 19.3 t ha −1 .

Simulation of silage maize
The observed LAI NDVI for silage maize in soil units A1ad (Figure 5a) showed rather homogeneous values on 28 May (LAI ∼ 0.4) and on 9 June (LAI ∼ 1.6). In contrast, the observed LAI NDVI for soil units D1a-d (Figure 5c) was more variable on 28 May (ranging from 0.3 to 0.8). This suggested a spatially heterogeneous presence of additional plants such as weeds within soil units D1a-d, which are not

F I G U R E 5 Simulations for silage maize. Simulated leaf area index (LAI, lines) and observed LAI calculated from normalized difference vegetation index (LAI NDVI , points) (a) in soil units A1a-d, (b) in soil units C1a, and (c) in soil units D1a-d (confidence interval of LAI NDVI shown with whiskers bars). Simulated water stress factor (d) in soil units A1a-d, and (f) in soil units D1a-d.
Note that a water stress factor of 1 does not produce a reduction of biomass accumulation (no water stress), whereas a factor of 0 represents the strongest simulated water stress. (e) Simulated mass of above surface organs in soil units A1a-d and D1a. Letters along the x axis represent months of the year, from March to November. Non-str., nonstressed simulated in AgroC. Afterwards, the observed LAI NDVI on 9 June (LAI ∼ 0.6) was lower than in May and also lower than that of the soil units A1a-d. This suggested the presence in units D1a-d of juvenile silage maize crops that were not present in May and that are in an earlier growth stage than those of soil units A1a-d. Therefore, simulations were performed using the same crop parameterization but with different emergence dates. The emergence date was set to 1 May for soil units A1a-d and to 10 May for soil units D1a-d.
The simulated LAI for soil units A1a-d ( Figure 5a) suggested homogeneous growth until late May. Afterwards, the simulated LAI of these four soil units diverged and showed considerable differences from June to August. The simulated LAI matched the observed LAI NDVI well, in general, although the simulated LAI of soil units A1b-c somewhat underestimated LAI NDVI . Since water stress was relatively low in June compared with the second part of the growing season, this variability in simulated LAI was attributed to the combination of the early stress period in May with later stress in July and August. The following abrupt decrease in LAI NDVI was due to the senescence stage of maize development. Interestingly, differences in LAI NDVI due to water stress were still present at this development stage. In the AgroC model, senescence in silage maize starts when a certain value of crop development stage is reached. In order to match the observed variation in LAI NDVI on 8 September, four different development stages were used to start senescence at four different times in the soil units A1a-d. Earlier senescence was assumed in soil units where stronger water stress was observed, since it was apparent from field observations how stronger stress resulted in smaller crops and early senescence. This required a modification of the model code that is described in Appendix C.
The simulated LAI for soil units D1a-c (Figure 5c) showed similar development throughout the growing season. This is consistent with the observed values of LAI NDVI . In contrast with soil units A1a-d, the water stress simulated in May was rather low for these three soil units (see water stress factor in Figure 5f) and did not affect the simulated LAI in the following months. Possibly, the later emergence date for soil units D1a-c prevented water stress from strongly affecting LAI in the early stage of maize growth. To corroborate this hypothesis, we performed the same simulation on soil unit C1a that typically does not show water stress. The simulated LAI of silage maize grown on this soil unit is shown in Figure 5b and matched well with the simulated LAI of soil units D1a-c. This confirms that the water stress factor in the early growth stage of silage maize has a strong influence on simulated LAI. In the case of soil unit D1d, simulated LAI was strongly reduced by water stress and did not match LAI NDVI values. Again, this is was attributed to the very small area of this soil-crop combination (0.05 ha). Figure 5e shows the simulated aboveground biomass of silage maize in soil units A1a-d and C1a, which is used as a measure of agricultural yield, since silage maize is F I G U R E 6 Simulations for potato. Simulated leaf area index (LAI, lines) and observed LAI calculated from normalized difference vegetation index (LAI NDVI , points) (a) in soil units A1a-d (no data since no potatoes were cropped), (b) in soil units BC, and (c) in soil units D1a-d (confidence interval of LAI NDVI shown with whiskers bars). Simulated water stress factor (d) in soil units A1a-d (no data since no potatoes were cropped), and (f) in soil units D1a-d. Note that a water stress factor of 1 does not produce a reduction of biomass accumulation (no water stress), whereas a factor of 0 represents the strongest simulated water stress. (e) Simulated mass of above surface organs in soil unit BC and D1a-d. Letters along the x axis represent months of the year, from March to November. Non-str., nonstressed used for livestock feeding in this area. The simulated yield ranged from 12.5 to 17.1 t ha −1 in soil units A1a-d and was 18.7 t ha −1 in soil unit C1a. These simulated yield values are in good agreement with the results of Žydelis et al. (2018), who used a comparable methodology and found similar measured and simulated yield: ∼18.7 t ha −1 in healthy silage maize and ∼14.8 t ha −1 in stressed maize. It has to be noted that the lower observed and simulated yield in Žydelis et al. (2018) was caused not only by water stress, but also by lower average air temperatures.

Simulation of potato
Figures 6b-6c show the simulated LAI and the observed LAI NDVI of potato in soil units C1a-b, C2a-b, and D1a-d. The observed LAI NDVI only showed small variations, except on 12 August where moderate differences between soil units C1a-b or C2a-b (from ∼4.4 to ∼4.7) and soil units D1a-d (from ∼3.8 to ∼4.3) were observed. Overall, the AgroC model was able to capture LAI development of potato fairly well. However, Figure 6f shows that the simulated water stress factor in soil units D1a-d was rather strong and affected the simulated LAI (Figure 6c) and the simulated dry aboveground biomass (Figure 6e), to some extent. However, no clear correlation between variations in simulated LAI and LAI NDVI due to water stress was observed. This is attributed to irrigation of the potato field cropped in 2016. Regrettably, detailed information on the amount, timing, and location of this irrigation was not available. Therefore, it was not possible to meaningfully implement irrigation in our simulations. It is expected that irrigation increased the LAI NDVI in soil units where strong water stress was simulated with AgroC, thus preventing a clear expression of water stress factor in LAI NDVI .

Simulation of winter wheat
The AgroC simulations for winter wheat are shown in Figure 7. Similar to sugar beet, three different crop parameterizations were used. In order to capture observed difference in maximum LAI NDVI , it was necessary to calibrate the partitioning of the mass allocated to the stem and to the leaf as described in detail in Appendix C. As in the case of sugar beet, one parameterization was used for each group of soil units shown in Figures 7a-7c. Thus, the variability in simulated LAI between soil units A1a-d and between soil units D1a-d can only be the results of different soil parameterizations. As shown in Figure 7, water stress was rather low throughout the growing season for winter wheat and was limited to three distinct periods: late February and March, May, and finally July. It is important to note that the simulated stress factor in February and March is caused by high water content (saturation). In soil units A1a-d (Figure 7a), simulated LAI is consistent with observed LAI NDVI . Both simulated and observed LAI showed little variation between soil units A1a-c, and lower values for soil unit A1d on 28 May and 9 June. In soil units BC (Figure 7b), the observed LAI NDVI was generally higher than in the other soil units. Here, the LAI NDVI values between all soil units were rather similar throughout the year, and the simulated LAI matched the observed LAI NDVI well. A similar result was obtained for soil units D1a-d, although the variability between soil units was somewhat higher than in the BC soil units. In the simulations, the variability between soil units D1a-c and D2a was not clearly reproduced (Figure 7c) because of low simulated water stress. On the contrary, soil unit D1d showed lower simulated LAI and observed LAI NDVI on 28 May and 9 June compared with the other soil units.
Despite the lower magnitude compared with summer crops, water stress affected simulated dry mass of wheat grains at harvest for the different soil units (Figure 7e). This is due to the fact that the storage organs grow closer to the end of the growing season (June and July) and are thus more influenced by water stress than by LAI (Steduto, Hsaio, Fereres, & Raes, 2012). Here, nonstressed soil units produced 7.7 t ha −1 of dry grains, whereas the most stressed soil unit produced only 6.6 t ha −1 . These values are comparable with actual yields observed in other studies (Han & Yan, 2017;Káš, Mühlbachova, & Kusá, 2018;Liu, Ren, Gao, Yan, & Li, 2017) for winter wheat growing under no-and low-water-stress conditions.

Simulation of winter barley
The simulated and observed LAI for winter barley are shown in Figures 8a-8c. Here, the same crop parameterization was used for all simulations except for the emergence date, which was set to 10 December for soil units A1a-d and to 1 December for the BC soil units and for soil units D1a-d. These different emergence dates resulted in lower simulated LAI for soil units A1a-d compared with D1a-d on 20 March, which is consistent with the observed LAI NDVI . Throughout the rest of the growing season, the variability in LAI NDVI between the soil units was rather strong but obviously not related to simulated water stress in soil units  or to the specific characteristics of each soil unit in BC. In summary, the simulated LAI of winter barley did not match the observed LAI NDVI as well as in the simulation of other crops but reproduced the general trend of LAI NDVI development after selecting appropriate emergence dates. One reason for this is that LAI NDVI showed differences between soil units, which are probably correlated to other factors than water stress such as crop management.
As discussed in the case of wheat, barley yield is more sensitive to changes in growing conditions during the period when the grain number is set (Steduto et al., 2012). For this reason, the water stress simulated in May had an influence on the final grain yield (Figure 8e). The simulated dry weight of barley grains in soil units D1a-d and BC in field F11 was 6.9 t ha −1 . This simulated value for field F11 was identical to the harvest reported by the field owner in 2016.

Simulation of winter rapeseed
Simulations for winter rapeseed were performed using a single crop parameterization, but the emergence date was set to 10 November for soil units A1a-d and to 1 November for soil units BC and D1a-d based on observed LAI NDVI . Figures 9a-9c shows the simulated LAI and the LAI NDVI of winter rapeseed, whereas Figures 9d and 9f show the water stress factor simulated in soil units A1a-d and D1a-d. A single peak of water stress was present in May, and no other stress was obtained throughout the rest of the simulation. In the absence of extended periods of water stress, only small differences in simulated LAI were observed between different soil units. In contrast, some differences in observed LAI NDVI was found between soil units A1a-d, especially for 14 March and 20 April (Figure 9a). It must be noted that the simulated LAI in May shows variations that are in agreement with the observed LAI NDVI of March and April. Nevertheless, these results suggest that the simulations underestimate the variability observed in the LAI NDVI . From 9 June, both simulated and observed LAI did not vary substantially between soil units A1a-d and showed a similar value (LAI ∼ 6.1).
The simulated LAI for soil units C1a-b, C2a-b, and D1a-d (Figures 9b-9c) did not match the observed LAI NDVI at all times. Obvious differences between simulated and observed LAI were present at 14 March and at 20 April. In March, the LAI NDVI of soil units D1a-d showed considerable differences with LAI NDVI values ranging from 3.8 for soil unit D1a to 4.9 for soil unit D1c. This variability was not captured by the LAI simulations. The reason for this difference in LAI NDVI may be associated with factors that are not considered in our modeling approach (e.g., fertilization or different management of separated fields). In April, the simulated LAI considerably overestimated the observed LAI NDVI . Here, the low observed LAI NDVI may have been influenced by the blooming of winter rapeseed that occurs in April or May, which reduces crop greenness and chlorophyll concentration. This has an influence on the NDVI estimated from satellite remote sensing and may have resulted in lower LAI NDVI values. It is possible that the LAI NDVI values of soil units A1a-d were not affected by this because the late emergence date may have delayed the blooming stage. Overall, it is important to note that the simulated water stress factor had a rather low impact on winter rapeseed growth. Figure 10 compares the simulated LAI of the 80 unique soil-crop combinations with the LAI NDVI derived from the six RapidEye observations. Generally, the spatial pattern in the observed LAI NDVI is well reproduced and the simulated emergence, growth, and senescence of each crop are well timed. However, a closer inspection allows to identify several discrepancies. For example, differences in simulated LAI and observed LAI NDVI of silage maize (fields F13b, F41, and F42 in Figure 10) and sugar beet (fields F01, F05, F12, F13a, F44, and F46-51 in Figure 10) can be observed in August and September (Figures 10e-10f). In this case, the general average of each soil-crop combination is well captured (Figures 4a-4c and 5a-5c). However, a lower spatial variability in simulated LAI was found compared with the observations. As mentioned above, the geophysics-based soil map successfully captured the main spatial variability in soil properties at the square kilometer scale, and thus the main patterns in soil-controlled water stress and the resulting differences in crop development. However, due to the concept of soil units, it is not expected that the geophysics-based soil map can reproduce pixel-scale variability of LAI NDVI that occurs within a single soil-crop unit. A lower spatial variability of simulated LAI compared with observed LAI NDVI at the early growth stage of winter rapeseed was also observed (Figure 10a). This might be caused by the absence of early water stress in the simulation of winter rapeseed (Figure 9a), but it may also be related to different management of adjacent fields that is not considered in the methodology. The effects of local management are considered to be of secondary importance at the square kilometer scale and cannot be meaningfully simulated with typically available information. Similar considerations apply to the case of within-field variability associated with soil compaction and tractors lines, which was evident in the fields cropped with potato (F40 in Figure 10) and locally in sugar beet (fields F01, F05, F12, F13a, F44, and F46-51 in Figure 10). Overall, the agroecosystem simulations for the growing season 2016 were capable of reproducing the growth of the six investigated crops, and the simulated LAI is generally consistent with the LAI NDVI observed from space.

Crop yield simulations at the square kilometer scale
The simulated yield at harvest for each soil-crop combination of sugar beet, silage maize, winter wheat, and winter barley is shown in Figure 11. In these maps, a value of 100% represents crops that have grown under optimal water supply (potential yield) for the entire growing season. Lower percentages are shown in the case of water-limited growth and correspond with a lower dry weight at harvest of storage organs (or aboveground biomass in the case of silage maize) relative to optimal water supply. All four crops achieved a 100% yield for the simulations of the BC soil units (except maize, which was not grown on units BC). For these soils, crops did not experience water stress, since the shallow ground water table provided sufficient water to guarantee optimal growth in 2016. On the contrary, a lower simulated yield was apparent for the soil units with underlying coarse layers. In this case, the yield of each crop showed a reduction that is largely proportional to the intensity of water stress. The lowest yield was commonly observed for soils with a shallow 2C layer with a large amount of gravels.
The simulated dry weight of storage organs in sugar beet was strongly influenced by soil type and ranged between 58 and 77% in soil units A1a-d, and between 62 and 92% in soil units D1a-d (Figure 11a). Silage maize simulated in soil units A1a-d ( Figure 11b) had a similar range but higher average yield compared with sugar beet (from 67 to 89%). Interestingly, the relative silage maize yield was much higher than that of sugar beet in soil units D1a-c and showed F I G U R E 10 Comparison between observed LAI calculated from normalized difference vegetation index (LAI NDVI , above) and simulated leaf area index (LAI, below) on (a) 14 Mar., (b) 20 Apr., (c) 28 May, (d) 9 June, (e) 12 Aug., and (f) 8 Sept. 2016. Field codes are colored according to the crop type (blue for sugar beet, purple for maize, brown for potato, and yellow for winter rapeseed). Fields are shown in gray when the simulated crop is not yet cropped or has already been harvested Vadose Zone Journal F I G U R E 11 Maps of the simulated agricultural yield in 2016 of (a) sugar beet, (b) silage maize, (c) winter wheat, and (d) winter barley much lower variability (from 99 to 100%). For the winter crops, simulated yield was still influenced by soil type, but the impact was generally lower, as evidenced by the higher average relative yield as compared with the summer crops. In particular, the simulated yield ranged from 86 to 94% for winter wheat (Figure 11c) and from 81 to 98% for winter barley (Figure 11d).
The maps presented in Figure 11 have valuable practical applications in the maximization of agricultural productivity. For example, the selection of a specific crop type could be based on the characteristics of the soils that are found in a given field. Furthermore, the seeding date of specific crops could be selected after investigating the precipitation recorded in the previous weeks and the forecasted meteorological conditions. Moreover, the proposed maps and simulation strategy could provide valuable information on the amount and timing of irrigation that is needed to reduce water stress and maximize crop yield.

SUMMARY AND CONCLUSIONS
In this study, the growth of six crops in 2016 within a 1-km × 1-km area was simulated using the agroecosystem model AgroC using data from a high-resolution geophysicsbased soil map. Soil hydraulic parameters were estimated from soil texture and bulk density using the PTF of Rawls and Brakensiek (1985). First, simulated water content of two soil units was compared with measured soil water content at three depths. A low RMSE of 0.032 cm 3 cm −3 was obtained when the soil profile provided by the geophysics-based soil map well represented the actual soil profile. Additional simulations performed using soil hydraulic parameters estimated using the ROSETTA PTF resulted in a much higher RMSE between measured and simulated water content, which was attributed to the higher estimated values of the saturated hydraulic conductivity. Future research could explore the use of more automated methodologies for the estimation of hydraulic parameters (e.g., inversion strategies), as well as the analysis of uncertainty in the simulation of soil water content dynamics.
In a next step, agroecosystem simulations were performed for six crops: sugar beet, silage maize, potato, winter wheat, winter barley, and winter rapeseed. In general, it was found that the magnitude of simulated water stress was a function of the crop type and of the soil characteristics. Higher water stress occurred in coarser soils and summer crops. Overall, the simulated LAI was found to be consistent with six maps of observed LAI NDVI , which were derived from RapidEye satellite data by using vegetation indices and ground truth information. Inconsistencies between simulations and observations were locally present (e.g., during the blooming of winter rapeseed, in the case of irrigation, or when the effect of field-scale management was not reproduced).
Water stress had an impact on the simulated yield at harvest for each of the investigated crops. Simulated yield of sugar beet and winter barley matched the actual harvest in 2016 in two fields within the study area, and the simulated yield of silage maize and winter wheat corresponded to literature values. Maps of the simulated yield were produced for these four crops, which showed an expected agricultural yield of 100% in soil units with a shallow ground water table. In soil units with underlying coarse soil texture layers, simulated yield showed a reduction that was largely proportional to the intensity of water stress. The simulated yield ranged between 59 and 92% for sugar beet, between 67 and 99% for silage maize, between 86 and 94% for winter wheat, and between 81 and 98% for winter barley.
Overall, this study showed that quantitative spatial information on soil heterogeneity derived from geophysics-based soil mapping allowed precise agroecosystem simulations of multiple crops on a large 1-km × 1-km area where water stress is strongly influenced by soil characteristics. Future studies could be focused on the quantification of the added value of geophysics-based soil mapping in comparison, for example, with commonly available and general-purpose soil maps. The thematic maps produced with the results of these simulations are relevant in practical agricultural applications, such as the selection of a specific crop type for a given soil or the selection of the best seeding date for crops that are particularly susceptible to water stress in the early growth stage. For this, additional information could be obtained by considering precipitation recorded in previous weeks or precipitation forecasted by weather models. By extending the simulation period (e.g., 30 yr), the proposed strategy could also allow the evaluation of the cost-benefit ratio of long-term strategies such as the most suitable crop rotation. Finally, the proposed maps and simulation strategy might provide valuable insights about irrigation scheduling and quantity by implementing information from forecasted precipitation.
How to cite this article: Brogi C, Huisman J, Herbst M, et al. Simulation of spatial variability in crop leaf area index and yield using agroecosystem modeling and geophysics-based quantitative soil information. Vadose Zone J. 2020;19:e20009. https://doi.org/10.1002/uar2.20009

APPENDIX A
The LAI NDVI maps used in this study were produced using Equation 4 and following the methodology described in Ali et al. (2015). For this, the NDVI was calculated for each image pixel using where NIR is the near-infrared (760-850 nm) and RED is the red spectral band (630-685 nm). In a next step, the FVC NDVI of each image was calculated using where NDVI s is the NDVI for bare soil and NDVI v is the value at the fully vegetated state (Beck, Atzberger, Høgda, Johansen, & Skidmore, 2006;Xiao & Moody, 2005;Zeng, Rao, DeFries, & Hansen, 2003). The values of NDVI s and NDVI v were estimated through histogram evaluation (NDVI s = −0.05 and NDVI v = 0.81).

APPENDIX B
In this appendix, a short description of the agroecosystem model AgroC is provided. This model couples (a) the SOILCO2 module (Šimůnek & Suarez, 1993;Šimůnek et al., 1996), (b) the RothC module (Coleman & Jenkinson, 1996), and (C) the SUCROS module (Spitters et al., 1989). The SOILCO2 module solves the one-dimensional Richard's equation that describes water flow in a given soil profile: where h (cm) is the pressure head, θ w (cm 3 cm −3 ) is the volumetric water content, K(h) (cm h −1 ) is the hydraulic conductivity as a function of pressure head, t is time (h), z is the vertical coordinate (cm), and Q (cm 3 cm −3 h −1 ) is the source-sink term accounting for water uptake by plant roots. The water retention and the unsaturated hydraulic properties as a function of pressure head are described by the Mualemvan Genuchten (van Genuchten, 1980) model: and where θ r and θ s (cm 3 cm −3 ) are the residual and saturated water content, respectively, K s is the saturated hydraulic conductivity (cm h −1 ), S e is the relative saturation (dimensionless), α is the inverse of the air entry pressure (cm −1 ), n is a dimensionless parameter related to the pore size distribution, and the parameter m is set equal to 1 − 1/n. In the SUCROS module, potential evapotranspiration of a crop (ET p,crop ) growing under optimal conditions is calculated by multiplying the potential grass reference evapotranspiration with a crop-specific coefficient K c . Here, the potential grass reference evapotranspiration is calculated from meteorological data using the Penman-Monteith approach (Allen et al., 1998). The value of the K c coefficient varies throughout the growing season with crop development stage. Then, ET p,crop is split into potential soil evaporation E p (cm h −1 ) and potential transpiration T p (cm h −1 ) according to Beer's law: where LAI is the green leaf area index and E i (cm h −1 ) is the water removed by evaporation of intercepted rainfall. The latter is calculated using where C i (cm) is the interception storage at a specific time step and S i (cm) is the canopy interception storage capacity that is assumed to be proportional to the LAI. The potential root water uptake S p (cm 3 cm −3 h −1 ) is a function of the potential transpiration T p over depth z according to the relative root density distribution. Actual root water uptake is calculated from the potential root water uptake using where φ (-) is a root water uptake stress factor that reduces the potential root water uptake in dry and wet soil conditions. This reduction factor is calculated using Equation 5 (Feddes et al., 1978). Finally, the actual transpiration is provided by the integration of the actual root water uptake over depth: where L r is the rooting depth (cm). At crop emergence, an initial rooting depth is specified for each crop. Thereafter, rooting depth increases in dependence of solar radiation, temperature, biomass production, and partitioning to the rooting system until the maximum rooting depth is reached. The root length density distribution with depth was calculated from a dimensionless weighting factor specified for relative rooting depths (Schamschula et al., 1974). Biomass production is estimated from the C assimilated by the plant, which mainly depends on temperature and solar irradiation. Here, the potential assimilation is scaled by the ratio T a /T p to also account for water stress. The glucose assimilation rate A g (kg CH 2 O m −2 h −1 ) is equal to about 30/44 of the total CO 2 assimilation rate A (kg CO 2 m −2 h −1 ), which is obtained by integrating the total instantaneous assimilation rate A L,T (kg CO 2 m −2 leaf surface h −1 ) over the leaf area. In a second step, the net growth rate of DM per unit area (GWT, kg DM m −2 h −1 ) is calculated using where ASRQ (kg DM kg −1 h −1 ) is a coefficient that represents the organic-specific conversion efficiency from glucose to DM. The R m (kg CH 2 O m −2 h −1 ) is the total maintenance respiration demand, which is controlled by plant senescence and temperature. After determining the net growth, the produced biomass is partitioned to the crop organs to describe DM accumulation in the roots, storage organs, stems, and leaves. This partitioning is crop specific and a function of development stage. During the crop juvenile stage, the development of LAI is driven by temperature. At later stages, the DM growth rate of green leaves (GWLV G , kg DM m −2 soil h −1 ) at a specific time step (t) is limited by the supply of assimilates. Therefore, the green LAI growth rate (GLAI G , m 2 leaf m −2 soil h −1 ) is where SLA (m 2 leaf kg −1 DM) is the specific leaf area. The leaf area growth rate decreases during the season because of senescence and self-shading. In the model, a dead LAI growth rate is calculated for different crops and is then subtracted from the GLAI G to describe this reduction in the growth rate.

APPENDIX C
In this appendix, additional information regarding the parameterization of sugar beet, silage maize, and winter wheat is provided. In sugar beet, a strong reduction of LAI NDVI in soil units A1a-d and in D1a-d was observed in August and September 2016 (Figure 4). At this stage, the AgroC simulations for LAI already reached their highest possible value, which was similar in all soil units. Here, the water stress simulated in May was not able to influence simulated LAI, since temperature is the main driver of LAI accumulation in the juvenile stage. At the same time, there was no possibility to reduce LAI in August and September by water stress because the simulated water stress does not feed back to the simulation of the death rate of leaves or senescence in the current AgroC model. However, the observed LAIs indicated the existence of this feedback mechanism. Therefore, it was decided to reproduce the LAI reduction in the late growing season by using three specific death rates of leaves (m 2 leaf m −2 soil • C d −1 ). The adopted death rate of leaves was generally higher in soil units A1a-d, intermediate in D1a-d, and lowest in D1ad (Table A1). This calibration influenced the simulated LAI in a way that it improved the simulation of the differences in observed LAI between subareas A, BC, and D. Since all other plant parameters and boundary conditions were identical in all simulations for sugar beet, the differences in simulated LAI within each subarea are solely due to the magnitude and timing of simulated water stress. In silage maize, the abrupt decrease in LAI that occurred in soil units A1a-d in September is due to the start of the senescence stage. Since the senescence stage can be affected by water stress (Baret et al., 2007), we assumed that crops subjected to greater water stress showed earlier senescence. This assumption was corroborated by field observations. In the AgroC model, senescence of silage maize is generally started when a development stage value of 1.4 is reached. To implement variable senescence consistent with stress intensity, this development stage value was set to 1.38, 1.35, 1.33, and 1.30 for soil units A1a-d, respectively.
In winter wheat, the three groups of soil units A1a-d, BC, and D1a-d showed differences in LAI NDVI , especially in April. In this case, the partitioning of aboveground biomass allocated to the stem and to the leaves was calibrated, and one specific partitioning was used for each group of soil units (see Table A2). Since one parameterization was used within each group of soil units, the variability in simulated LAI in each group is the results of different soil parameterizations that affect water stress.