Can big data explain yield variability and water productivity in intensive cropping systems?

Yield gaps and water productivity are key indicators to monitor the progress towards more sustainable and productive cropping systems. Individual farmers are collecting increasing amounts of data (‘big data’), which can help monitor the process of sustainable intensification at local level. In this study, we build upon such data to quantify the magnitude and identify the biophysical and management determinants of on-farm yield gaps and water productivity for the main arable crops cultivated in the Netherlands. The analysis focused on ware, seed and starch potatoes, sugar beet, spring onion, winter wheat and spring barley and covered the period 2015–2017. A crop modelling approach based on crop coefficients (kc) and daily weather data was used to estimate the potential yield (Yp), radiation intercepted and potential evapotranspiration (ETP) for each crop. Yield gaps were estimated to be ca. 10% of Yp for sugar beet, 25–30% of Yp for ware, seed and starch potato and spring barley, and 35–40% of Yp for spring onion and winter wheat. Variation in actual yields was associated with water availability in key periods of the growing season as well as with sowing and harvest dates. However, the R2 of the fitted regressions was rather low (20–49%). Current levels of crop water productivity ranged between 13 kgDMhamm for spring barley, ca. 15 kg DMha−1mm−1 for seed potato, spring onion and winter wheat, 23 kg DMha−1mm−1 for ware potato and ca. 25 kg DMha−1mm−1 for starch potato and sugar beet. These values are about half of their potential, but increasing actual water productivity further is restricted by rainfall amount and distribution. However, doing so should not be prioritized over reducing environmental impacts of these intensive cropping systems in the short-term and may require large investments from farm to regional levels in the long-term. Although these findings are most relevant to similar cropping systems in NW Europe, the underlying methods are generic and can be used to benchmark crop performance in other cropping systems. Based on this work, we argue that ‘big data’ are currently most useful to describe cropping systems at regional scale and derive benchmarks of farm performance but not as much to predict and explain crop yield variability in time and space.


Introduction
Sustaining food production in the years ahead requires sustainable and efficient management of natural resources and the surrounding environment. Different indicators have been suggested to measure the progress towards more sustainable and productive agricultural systems (Gil et al., 2019). The yield gap is an important indicator in this context as it provides a benchmark for agricultural productivity but its interpretation must go alongside other indicators expressing the efficiency of critical resources such as water and nutrients. Benchmarks for the latter include water productivity or nitrogen (N) use efficiency and N surplus.
weeds. The aforementioned yield levels reflect differences in crop performance in terms of resource use efficiency. For instance, Yw provides a benchmark for water productivity (WP) which is defined as the dry matter produced per mm of water available.
Both yield gap and water productivity are important indicators to benchmark agricultural systems . Traditionally, these have been assessed with crop growth models considering the variation in biophysical conditions and, to a lesser extent, crop management practices observed in farmers' fields. The advent of individual farmer field data brings new opportunities to assess crop performance under farmers' conditions and to measure the on-farm impact of technological innovations over large geographic areas (e.g., Berrueta et al., 2019;Rattalino-Edreira et al., 2018;Silva et al., 2017b;Delmotte et al., 2011). In fact, they 'can be considered equivalent to running hundreds of field experiments to capture both major management effects and management x environment interactions' (Rattalino-Edreira et al., 2018). Its use involves skills and methods for the manipulation of 'big data' (i.e., high volume, velocity and variety of information to require specific analytical and technological methods for its transformation into value) beyond those used in traditional agronomic research (Wolfert et al., 2017;de Mauro et al., 2016).
Farm level data has been collected in the European Union for a long time (e.g., Farm Accountancy Data Network) to monitor and evaluate the income of agricultural holdings and changes in farm structure. These were also used to investigate adaptation options to climate change at farm level in the European Union (Reidsma et al., 2010) or to decompose yield gaps of the main arable crops in the Netherlands (Silva et al., 2017a). However, few spatially explicit databases exist with detailed biophysical and crop management information for a large number of farmers' fields (Beza et al., 2017). Individual farmers are recording an increasing amount of data but it remains unclear how these can help gaining deeper insights into the drivers of resource use efficiencies and accelerating the process of sustainable intensification at local level.
The objectives of this study are to quantify the magnitude and identify the drivers of yield gaps and water productivity of key arable crops cultivated in Northwest Europe, using the Netherlands as a case study. We hypothesize yield gaps and water limitations are small in farmers' fields but resource use efficiencies can be improved (cf. Silva et al., 2017a). The analysis draws upon a large database of individual farmer field data during the period 2015-2017 and covers the production ecology of root, tuber and bulb crops (ware, seed and starch potato, sugar beet and onion) as well as cereal crops (winter wheat and spring barley). The benchmarks and insights presented in this study help understanding the scope to improve productivity and resource use efficiency of intensive arable cropping systems in the Netherlands and other similarly intensive cropping systems.

Dutch arable farming
Arable farms in the Netherlands cultivate different types of crops, the most important ones being ware, seed and starch potatoes, sugar beet, spring onion, winter wheat and spring barley (CBS, 2015). These crops are often integrated into crop rotations in which each crop is recurring in the same field every certain number of years. Yield gaps in these cropping systems are amongst the smallest worldwide, with Ya reaching 70-80% of Yp Silva et al., 2017a;Rijk et al., 2013). These small yield gaps are the result of modern varieties and cultivation techniques as well as of intensive use of inputs, particularly fertilizers (mineral and organic) and pesticides. Irrigation is not a default in the country, due to the humid climate and shallow depth of groundwater, but it is getting increasing attention as a climate adaptive strategy to dry summers (e.g., Kroes et al., 2018;van Duinen et al., 2015). This justifies using Yp as a benchmark instead of Yw (Silva et al., 2017a).
The Netherlands has a flat topography and three dominant soil types for arable cropping: marine clay soils in the north and south-western areas, sandy soils in south-and central-eastern areas and, less representative, peat soils in the north and center-western areas. Climatic conditions are rather homogeneous throughout the country, which is characterized by a temperate maritime climate (Fig. 1). Solar radiation, minimum and maximum temperatures and reference evapotranspiration increase in the first half of the year to a maximum value per month of ca. 17.5 MJ m −2 , 12.5°C, 22.5°C and 110 mm, respectively, during the summer, after which there is a decline in the second half of the year. Rainfall is ca. 800 mm yr −1 and it is distributed fairly evenly throughout the year.

Background and spatial distribution
The database used in this study included more than 5600 field × year combinations. Data were self-recorded by farmers in farm management systems of a consultancy company providing crop registration software and decision support systems in the Netherlands. Farmers use these software platforms to register biophysical characteristics and crop management operations at field level including sowing and harvest dates, water and nutrient management, pest and disease management and final crop yield.
The relational database with the raw data contains one table per management operation, which can be linked through an unique farm × field identifier. The management operations recorded are land preparation (date and method), crop establishment (sowing date, density, method and variety), irrigation (dosage, date, source and method), fertilization (idem), pesticide use (idem) and harvest (crop yield, date and method). Two additional tables include farm-and field-specific information such as location (postal code and GPS coordinates), commercial software used, crop year (2015-2017), cultivated crop, field area (ha) and soil type (based on the Dutch soil classification; cf. Hartemink and Sonneveld, 2013). Soil physical (e.g., texture) and chemical properties (e.g., organic matter, pH and available NPK) from soil samples are also available for a smaller number of fields (ca. 20% of all fields).
The fields included in the database are located across the main agricultural regions of the Netherlands (data not shown due to privacy reasons). Cultivation of potato is spread across the country. Ware potato for the French fries industry and direct consumption is mostly produced in the polders (i.e., Flevoland) and southern regions, seed potato is mostly found in the northeast and north coastal regions and starch potato is confined to the northeastern region. Sugar beet cultivation is spread throughout the country, while spring onion is mostly produced in the southwestern, polders and northern coastal regions. Winter wheat is cultivated in similar regions as spring onion, while spring barley is mostly found in the northeastern region. Wheat is mainly produced for animal feed rather than bread baking while barley is produced for brewing beer and for animal feed.

Data management and secondary sources
Different queries were designed with the Python module pyodbc to retrieve from the database all data related to each management operation, with the exception of the pesticide use which was beyond the scope of this study. These queries were screened for basic problems (e.g., duplicated records and non-standard units) and used to construct the variables of interest (i.e., total input applied and input dose, date, method and product per split). Actual yields (Ya) were retrieved from the database and standardized to ton fresh matter per ha (t FM ha −1 ). Soil types classified as 'loess', 'Maas clay', 'old marine clay', 'peat', 'peat (with clay)', 'river clay' and 'young marine clay' in the Dutch classification system were grouped as clay soils (n > 3200) and soil types classified as 'cover sand', 'reclaimed peat soil' and 'marine sand' were grouped as sandy soils (n > 2400). Potato varieties were classified based on maturity ('early', 'medium' and 'late'), sugar beet varieties based on disease/nematode resistance ('not resistant' and 'resistant' against rhizoctonia, rhizomania, cyst nematodes), onion varieties based on color ('yellow' and 'red'), wheat varieties as 'bread-baking' quality or 'other purposes ' and, barley varieties as 'malt' or 'feed' barley. Daily weather data were obtained from the Royal Dutch Meteorological Agency (KNMI, in Dutch) for 25 weather stations spread across the Netherlands. These contained complete records of global solar radiation (MJ m −2 ), minimum and maximum temperature (°C), precipitation (mm), vapour pressure (hPa) and wind speed (m s −2 ) during the period of interest (2015-2017, Fig. 1). Reference evapotranspiration (ET0, mm) was estimated for each weather station using the FAO-Penman-Monteith equation (Allen et al., 1998). The nearest weather station to each farm was identified based on the postal code of the farm as GPS coordinates were not available for all fields.

Crop modelling approach
A framework based on crop coefficients (kc) was used to estimate the potential yield (Yp), radiation intercepted and potential evapotranspiration of ware, seed and starch potato, sugar beet, spring onion, winter wheat and spring barley (Fig. 2). This summary model was preferred over a mechanistic crop growth model due to the large Crop modelling approach used to simulate Yp, crop water requirements and soil water losses for arable crops in the Netherlands. The crop coefficients (kc**) used were adjusted to the length of the growing season of each field -see text for further explanation. Rate equations are provided as Supplementary Material and the parameters used for each crop are provided in Table 1. See also Sadras et al. (2017) for further details. Abbreviations: P% = protein content, F% = fat content, HI = harvest index, RUEp = potential radiation use efficiency, RUEc = reference RUE, f_PAR = fraction of intercepted PAR. number of crops analysed and the lack of up-to-date parameters for some arable crops in the Netherlands to run the latter. Soil water losses and seasonal water available were also estimated using a daily water balance based on a 'tipping bucket' approach for two contrasting soil types (clay and sand). A detailed description of the approach is provided below.

Yield potential and crop water requirements
The crop modelling approach used to estimate Yp of each crop is summarized in Sadras et al. (2017). This approach defines Yp as a function of solar radiation and crop characteristics and builds upon daily records of global radiation (RAD GLOBAL , MJ m −2 day −1 ) and a number of crop-specific parameters ( Table 1). The latter include harvest index (kg kg −1 ), dry matter (DM%), protein (P%) and fat (F%) contents (g protein or fat g DM −1 ) as well as tabulated crop coefficients (kc_ini, kc_mid and kc_fin, mm mm −1 ) and crop growth periods (A, B, C and D days).
The crop coefficient (kc) is defined as the ratio between the actual evapotranspiration of a given crop (ETc, mm) and ET0. The value of kc during the growing season was estimated according to Doorenbos and Pruitt (1977), which represents the kc curve as a set of straight lines. This curve can be defined based on the duration of different crop development stages (A, B, C and D), which in total define the length of the growing season (length A+B+C+D ), and the value of kc at three points in time (kc_ini, kc_mid and kc_fin). Stage A refers to the initial period between sowing and ca. 20% ground cover, stage B refers to the period of rapid crop growth from 20% to 70-80% ground cover, stage C refers to the period of maximum and constant ground cover and, finally, stage D refers to the period in which ground cover declines due to crop senescence until the crop is harvested or reaches physiological maturity . The value of kc was assumed to be constant during phases A and C and equal to kc_ini and kc_mid, respectively. During phases B and D, the value of kc was estimated through linear interpolation between kc_ini and kc_mid and between kc_mid and kc_fin, respectively. This approach assumes an average growing season, which is inadequate to estimate Yp for a large set of farmers' fields. To overcome this limitation, correction factors were used to increase or decrease the average duration of each phase (Av-Length, Table 1), depending on the length of the growing season of each field (Length FIELD ) as per Equation A1 (see Supplementary Material). These correction factors were calculated as the difference between the maximum and minimum duration for each crop stage, divided by the difference between the maximum and minimum length of the growing season (i.e., the sum of maximum or minimum durations of A, B, C and D). Maximum and minimum durations of crop stages were obtained from Villalobos and Fereres (2017), who reported values adapted from Doorenbos and Pruitt (1977) and Allen et al. (1998).
Total aboveground biomass production (TAGP, g DM m −2 , also including below-ground storage organs in case of root, tuber and bulb crops) under potential conditions was defined as a linear function of the PAR intercepted during the growing season (PAR INT , MJ PAR m −2 ), considering a crop-specific radiation use efficiency coefficient (RUEp, g DM MJ PAR m −1 ; Eq. A2; Sadras et al., 2017;Steduto et al., 2009;Monteith et al., 1977). PAR INT was calculated following Eq. A3 as the product of the daily global radiation (RAD GLOBAL ) and the fraction of intercepted PAR by green leaves on a given day (f_PAR), integrated between the sowing and haulm killing date for potato or harvest date for the other crops. Haulm killing of potato was assumed to have occurred three weeks before the reported harvest dates, which represents average farm practices in the Netherlands. The crop development stage is reflected in f_PAR by defining it as a function of kc (Eq. A4; Allen et al., 1998;Doorenbos and Pruitt, 1977). Incident PAR was assumed to be 0.45 of RAD GLOBAL (Sinclair and Muchow, 1999). Crop-specific RUEp for dry matter production was defined in relation to a reference value RUEc (assumed as 3.2 g DM MJ PAR −1 for C3 crops; Sadras et al., 2017) considering the crop-specific harvest index (HI), protein (P%) and fat contents (F%; Equation A5). These parameters are provided in Table 1 and were obtained from Sadras et al. (2017) and other secondary data sources. Yp (t FM ha −1 ) was calculated based on TAGP, HI and DM% (Eq. A6) and the yield gap closure was calculated as the ratio between Ya and Yp. The actual radiation use efficiency (RUEa) was further estimated as the ratio between Ya and PAR INT (Eq. A7). The aforementioned calculations were done for unique combinations of farm × field × year × crop in the database, using the farmer reported crop calendar and the global radiation data from the nearest weather station.
Crop water requirements during the growing season were estimated following the methodology described by Steduto et al. (2009). It was assumed that crop water requirements for potential production are equal to the cumulative potential crop evapotranspiration (PCETP, mm) during the growing season, which was calculated based on the daily kc and daily ET0 (Eq. A8). By definition, PCETP equals the sum of potential crop transpiration (TRP, mm) and total evaporation from soil and canopy (SEV, mm). The former depends on crop development Table 1 Crop-specific parameters used to simulate Yp, radiation intercepted and crop water requirements. Similar parameters were used for ware, seed and starch potato unless indicated differently. Sources of these values are provided in the main text. J.V. Silva, et al. Field Crops Research 255 (2020) 107828 (expressed by kc), weather conditions (expressed by ET0) and intercepted radiation (expressed by f_PAR; Eq. A9), while the latter was computed as the difference between PCETP and TRP (Eq. A10). PCETP comprises the amount of water needed to achieve Yp for each unique combination of farm × field × year × crop included in the database. For each crop, Yp and PCETP were used as benchmarks for Ya and seasonal water available, respectively, and to derive the potential water productivity (i.e., ratio between Yp and PCETP).

Seasonal water available and soil water losses
The total seasonal water available for crop growth (TSWA, mm), and actual water productivity (WP, kg DM ha −1 mm −1 ), were calculated for each unique farm × field × year × crop combination included in the database as follows: where θ sowing (mm) is the soil water available at sowing, RAIN (mm) is the effective growing season rainfall from the nearest weather station, IRRIG (mm) is the total irrigation water supplied, CapR (mm) is the amount of water available through capillary rise and DP (mm) is the amount of water lost through deep percolation (drainage). θ sowing was determined for clay and sandy soils as the difference between field capacity (FC) and permanent wilting point (PWP) times the rooting depth of the crop (Z; 0.5 m for root, tuber and bulb crops and 1.2 m for cereals), assuming soils were at 50% of the field capacity at sowing. Effective rainfall was equal to the rainfall between sowing and harvest for all crops except winter wheat, for which only rainfall after March 1st was considered. Average values of capillary rise were assumed per crop and soil type based on data from Kroes et al. (2018) and losses through deep percolation were estimated for each soil type following a daily soil water balance as described below. The outflows of water from the soil profile considered in the water balance were crop transpiration and losses through soil evaporation, deep percolation and surface run-off (Eq. A12). Soil water losses are most evident at the beginning of the growing season, when evapotranspiration is lower than rainfall (promoting drainage) and when soil coverage is insufficient to avoid direct evaporation losses. Soil water losses were estimated following a daily soil water balance that considered inflows of rainfall and capillary rise (Eq. A13). The water balance was simulated for a clay (FC = 0.27, PWP = 0.12 and saturation as SAT = 0.37, all in mm mm −1 ) and sandy soil (FC = 0.21, PWP = 0.10 and SAT = 0.31) using Eqs. A11-A13. Surface run-off was assumed to be negligible due to flat topography in the Netherlands. Evaporation was estimated with Eqs. A8-A10 and deep percolation was calculated for each soil type using a 'tipping bucket' approach (Eq. A14). In general, soil water losses were largely attributed to soil evaporation rather than to deep percolation (Fig. A4A).

Statistical analysis
Maximum boundary functions were estimated using quantile regressions fitted to the 90th percentile of the data. All quantile regressions were estimated per crop using the rq() or nlrq() functions of the quantreg package in R (Koenker, 2016). Different functional forms were specified for different independent variables. A quadratic functional form (y = ax 2 + bx + c) was assumed for the relationship between Yp or Ya, on the one hand, and sowing or harvest date, on the other. If the quadratic term was not statistically significant at 5%, then a linear functional form (y = ax + b) was fitted to the data. If the slope of the linear regression was also not statistically significant at 5%, then no quantile regression was fitted to the data. Conversely, a 'broken-stick' functional form (y = ax + b, x < c and y = d, x ≥ c) was estimated for the relationship between Ya and radiation intercepted, and between Yp and potential crop evapotranspiration (Gilbert and Hernandez, 2019). Finally, a logistic functional form (y = x + 0.99 x ) was estimated for the relationship between the radiation intercepted and the length of the growing season.
Multiple linear regression was used to assess the relationship between Ya and a set of biophysical and management factors. The biophysical variables considered were year, soil type and region dummies, intercepted PAR, and spring and summer rainfall. Management variables included sowing and harvest dates (linear and quadratic effects), variety type expressed as a dummy, spring and summer irrigation, total amount of plant available N, P and K applied and, field size. All variables included in the models were not correlated and had a variance inflation factor smaller than 10 (data not shown). The intercept and coefficients of the regression models were estimated for the pooled sample of each crop using ordinary least squares (OLS, lm function in R) with all continuous variables mean-scaled. The relationship between the fitted values and the residuals was checked for normality and homoscedasticity and no patterns were observed in the residual plots of the different models (data not shown).

Yield gaps in farmers' fields
Yield gaps were smallest for sugar beet (12% of Yp), intermediate for ware, seed and starch potato and spring barley (25-30% of Yp) and largest for spring onion and winter wheat (35-40% of Yp; Fig. 3A). For sugar beet, Ya and Yp were on average 86.5 and 99.0 t FM ha −1 . Consistent differences in Ya and Yp were observed for potato production systems: yields were greatest for ware potato (Ya = 52.7 and Yp = 72.2 t FM ha −1 ), intermediate for starch potato (Ya = 44.8 and Yp = 61.1 t FM ha −1 ) and smallest for seed potato (Ya = 36.6 and Yp = 54.2 t FM ha −1 ). For cereals, Ya and Yp were on average 9.6 and 15.3 t FM ha −1 for winter wheat and 6.7 and 9.8 t FM ha −1 for spring barley, respectively. The large yield gaps for spring onion were a result of low Ya (58.9 t FM ha −1 ) compared to Yp (90.5 t FM ha −1 ).
Highest farmers' yields (Y HF ) were estimated as the average top 10th percentile of Ya for unique crop × year × soil type combinations. These were 67.8, 46.3 and 52.8 t FM ha −1 for ware, seed and starch potato, respectively. Y HF of sugar beet was 103.6 and of spring onion 81.6 t FM ha −1 . For cereals, Y HF ranged between 8.5 and 11.1 t FM ha −1 for spring barley and winter wheat, respectively. For all crops, Y HF was greater in clay soils than in sandy soils (Figs. A1A and A1E). Differences in biophysical conditions and crop management between highest-(Ya ≥ 90th percentile), average-(10th percentile ≤ Ya ≤ 90th percentile) and lowest-yielding fields (Ya ≤ 90th percentile) were generally similar to those obtained with the multiple regression analysis explained below (Fig. A1).

Variety types
Ya of ware and starch potato were significantly greater for medium maturity varieties than for early maturity varieties and there were no significant yield differences between early and late maturity varieties (Table 2). For sugar beet, significantly greater Ya was observed for resistant varieties than for varieties not resistant to diseases or nematodes and the greatest Ya was observed for varieties resistant to both Rhizomanie and cyst nematodes. Ya was significantly greater for yellow varieties than for red varieties of spring onion and there were no significant yield differences between varieties for cereals.
The different maturity types of ware potato were sown on average in mid-April but medium and late maturity types were harvested ca. 25 days later than early maturity types (around mid-October compared to mid-September; Table 3). This difference in growing season resulted in greater cumulative radiation intercepted (670 vs. 610 MJ PAR m −2 )  and rainfall (445 vs. 410 mm) for the medium and late maturity varieties, but in a similar radiation use efficiency (ca. 1.8 g DM MJ PAR −1 ) and water productivity (ca. 22 kg DM ha −1 mm −1 ), compared to early maturity varieties. The same was true for seed potato, although differences in harvest dates between maturity types were less pronounced than for ware potato as seed potato is harvested immature depending mostly on the quality class of the seed. Sugar beet varieties not resistant to diseases or nematodes were used only in few fields (n = 5) but these tended to have lower yield and lower water productivity than resistant varieties (Table 3). Varieties resistant to Rhizoctonia were harvested ca. 10 days later than other types of resistant varieties. Crop calendars and resource use efficiencies were similar for different varieties of sugar beet, spring onion, winter wheat and spring barley.

Biophysical conditions
Significant differences in Ya were observed across regions for all crops (Table 2). On average, Ya of sugar beet and winter wheat was greatest in Flevoland, Ya of ware potato was greatest in the north coast while for seed potato the greatest Ya was observed in the south, north west, east and north coast regions. Ya of spring onion was significantly lower in the south west compared to the polders, while Ya of spring barley was lower in the east, south and south west. There were also significant yield differences across years for all crops, except spring onion ( Table 2). Ya of ware and seed potato was significantly greater in 2017 and 2016 than in 2015, respectively. For starch potato and winter wheat, Ya was significantly smaller in 2016 than in 2015. Significantly greater Ya of sugar beet was observed in 2017 compared to 2015, while the opposite was true for spring barley. Ya of seed potato and sugar beet was 2-3 t FM ha −1 lower in sandy soils than in clay soils, while for winter wheat this difference was about 0.5 t FM ha −1 .
Summer rainfall had a significant positive effect on the yields of ware potato, spring onion and winter wheat and spring rainfall had a significant negative effect on sugar beet yields (Table 2). Supplementary irrigations during spring and summer were also positively associated with ware potato yields. Although the effects of rainfall and irrigation on crop yields were small, these results indicate that the time at which water is available to crops is as important as the amount. Intercepted PAR was positively associated with cereal yields, and negatively associated with ware potato yields, but in both cases the effects were small.

Crop management
The effects of sowing and harvest dates on crop yield were clearer for all crops, except spring onion, when looking at the highest-yielding fields for a given sowing date (Figs. 3, A2 and A3) rather than at average-yielding fields (Table 2). Simulated Yp in the year 2017 decreased linearly with delayed sowing for all crops, except winter wheat for which there was no significant effect of sowing date on Yp (Fig. 3). The effects of sowing dates on Ya observed in highest-yielding fields, also in 2017, were similar to those observed for the Yp. Delaying sowing by one day reduced Ya in highest-yielding fields by 340 kg FM ha −1 day −1 for seed potato, 910 kg FM ha −1 day −1 for sugar beet, 20 kg FM ha −1 day −1 for winter wheat and 80 kg FM ha −1 day −1 for spring barley (Fig. 3). The effect of later sowing on Ya was less clear for ware and starch potato than for the aforementioned crops. Finally, the sowing date of root, tuber and bulb crops was not associated with the length of the growing season, while for cereals later sowing resulted in a significantly shorter growing season (data not shown).
The harvest window of root and tuber crops was longer than for cereals (Fig. 3) and the effects of harvest date on Ya and Yp were consistent across years (Figs. 3, A2 and A3). During 2017, there was a significant quadratic relationship between harvest date and the maximum Yp simulated for a given harvest date for tuber, root and bulb crops (Fig. 3A-F). For winter wheat, the maximum Yp for a given harvest date increased linearly with harvest date (Fig. 3G) while there was no relationship between both for spring barley (Fig. 3H). The effect of harvest date on the Ya observed in highest-yielding fields was different for different crops. For instance, for ware potato and sugar beet there was a quadratic relationship with a maximum between harvest date and Ya in highest-yielding fields ( Fig. 3B and E) while for seed potato and winter wheat there was a negative linear relationship between both (Fig. 3C and G).
There was a positive significant effect of plant available N applied on Ya of ware and seed potato, winter wheat and spring barley, but again the effect was rather small (Table 2). Application rates of P were also positively associated with seed potato and winter wheat Ya and application rates of K were positively associated with seed potato and spring barley Ya. These effects were also rather small except for P in seed potato for which Ya increased by 160 kg FM ha −1 kg −1 P applied. Overall, these results suggest that yield responses to N, P and K were rather small. Finally, increasing field size by 1 ha resulted in an increase of ware potato and sugar beet Ya of ca. 300 kg FM ha −1 and starch    (Table 2). A small positive effect of field size on Ya was also observed for cereals.

Radiation intercepted by arable crops
Total radiation intercepted during the growing season ranged between 400 and 900 MJ PAR m −2 for root, tuber and bulb crops, 700 and 950 MJ PAR m −2 for winter wheat and, 500 and 700 MJ PAR m −2 for spring barley (Fig. 4A-C). Radiation intercepted was linearly related to crop yield up to 660-680 MJ PAR m −2 for spring barley, spring onion and potato crops, 740 MJ PAR m −2 for sugar beet and 880 MJ PAR m −2 for winter wheat. After these values, there were no further increases in Ya with increasing amount of radiation intercepted. The slope of the quantile regression ranged between 1.2 g DM MJ PAR −1 for cereals, 2.3 g DM MJ PAR −1 for potato and sugar beet and, 2.4 g DM MJ PAR −1 for spring onion. Radiation intercepted increased non-linearly with the length of the growing season especially for potato crops and sugar beet ( Fig. 4D-E). For spring onion, spring barley and winter wheat, radiation intercepted increased in a more linear fashion with increases in the length of the growing season (e.g., Fig. 4F).

Benchmarking water productivity
Actual water productivity was greatest for sugar beet and starch potato with an average of ca. 25 kg DM ha −1 mm −1 and lowest for spring barley with an average of 13 kg DM ha −1 mm −1 (Fig. 5A). Intermediate values of ca. 15 kg DM ha −1 mm −1 were estimated for seed potato, spring onion and winter wheat, and of 23 kg DM ha −1 mm −1 for ware potato. For all crops, these average values are much lower than the potential water productivity calculated as the ratio between Yp and PCETP (Fig. 5A): ca. 40 kg DM ha −1 mm −1 for ware potato, starch potato and sugar beet, ca. 35 kg DM ha −1 mm −1 for seed potato and spring onion, 28 kg DM ha −1 mm −1 for winter wheat and 22 kg DM ha −1 mm −1 for spring barley. The large water productivity gap resulting from the difference between these two metrics is a result of yield gaps and relatively large water surplus in most fields.
The x-intercept, slope, inflection point and plateau of the boundary functions of Yp as a function of PCETP are biophysically meaningful ( Fig. 5B-H; e.g., Lollato et al., 2017;Grassini et al., 2011;French and Schultz, 1984). The x-intercept indicates the soil water losses due to evaporation, and it was greater for sugar beet and spring onion (> 120 mm) than for the other crops. The slope of the regression is a measure of potential water productivity (expressed in kg yield per unit transpiration) and it had a value of 34 kg DM ha −1 mm −1 for spring barley, 42 kg DM ha −1 mm −1 for winter wheat, ca. 50 kg DM ha −1 mm −1 for ware, seed and starch potato, and ca. 63 kg DM ha −1 mm −1 for sugar beet and spring onion. The estimated inflection point indicates the minimum amount of water needed to reach Yp and it ranged between 310 mm for seed potato, 370 mm for spring barley, ca. 415 mm for ware potato, starch potato and winter wheat, and ca. 510 mm for sugar beet and spring onion. Beyond these values no water limitations are expected and hence, the plateau of the regression lines indicates the maximum Yp for each crop.
These biophysical benchmarks provide a reasonable upper limit for the water productivity estimated in farmer's fields based on Ya and seasonal water available (Fig. 5). Although the latter was greater than the inflection point described above for most crop × field × year combinations, meaning that overall there was enough water to reach Yp, we note that in the dry year 2016 (despite very high rainfall and flooded fields in June, Fig. 1B) nearly 50% of the ware potato fields, 75% of the starch potato fields and 60% of sugar beet fields may have experienced water stress due to insufficient water. We also note that fields where water was available in non-limiting amounts may also have experienced water limitations due to sub-optimal distribution of rainfall during the growing season (ca. 80% of seasonal water comes from rainfall, Fig. A4B). This may well explain why some fields achieved higher Ya than others with similar amounts of seasonal water available.
There was a significant positive linear relationship between the length of the growing season and the seasonal water available for every crop during most of the years, but the effects were less pronounced for winter wheat (Fig. 6). The slope of the linear regression was greatest in 2017 for seed potato (2.2 mm day −1 ), starch potato (1.5 mm day −1 ), spring onion (3.4 mm day −1 ) and winter wheat (0.7 mm day −1 ). These effects were more variable for the years 2015 and 2016. For instance, longer growing seasons did not translate into greater seasonal water available in 2015 for spring onion and in 2016 for seed potato while the slope of this relation was greatest in 2015 for ware potato (3.2 mm day −1 ). Although there were significant effects of growing season length on seasonal water available for winter wheat, these were small (0.5-0.7 mm day −1 ) in 2015 and 2017 and meaningless in 2016 (−0.4 mm day −1 ). For spring barley, the largest effect of growing season on seasonal water available occurred in 2016 (3.2 mm day −1 ), followed by 2015 (1.6 mm day −1 ) and 2017 (1.5 mm day −1 ). Finally, the length of the growing season explained more than 50% of the variation observed in seasonal water available for all crops (adjusted-R 2 not shown). The large variation in seasonal water available for a given length of the growing season is explained by differences in effective rainfall across the country (Table A1) and, to a lesser extent, by differences in soil water and capillary rise between clay and sandy soils.

Limited scope to narrow yield gaps further
Yield gaps in intensive arable cropping systems of Northwest Europe are among the smallest worldwide Mueller et al., 2012;Neumann et al., 2010). This is true for most arable crops (e.g., Schils et al., 2018;Beza et al., 2017;Silva et al., 2017a) due to greater yield progress in farmers' fields over the past decades than the genetic progress in Yp. Here, we estimated yield gaps in the Netherlands to be ca. 10% of Yp for sugar beet, 25-30% of Yp for ware, seed and starch potato and spring barley, and 35-40% of Yp for spring onion and winter wheat ( Figure 3A). Similar yield gaps for ware potato, spring onion and spring barley were reported in Silva et al. (2017a). For starch potato and sugar beet, yield gaps were ca. 15% smaller, and for winter wheat ca. 10% greater, than those estimated by Silva et al. (2017a). These differences can be explained by a greater Yp for starch potato, lower Ya for sugar beet and lower Yp for winter wheat in Silva et al. (2017b), which covered the period 2008-2012 and used a different crop growth model.
The Yp simulated in this study were generally greater than those reported by Silva et al. (2017a) for the same cropping systems. The Yp used by Silva et al. (2017a) were obtained from variety trials for starch potato, sugar beet and spring onion (Rijk et al., 2013) and from crop model simulations using WOFOST for ware potato, winter wheat and spring barley (Reidsma et al., 2015b). We note that these previous benchmarks were lower than the Ya in some of the fields analysed here, particularly for ware potato, sugar beet and winter wheat (Fig. 3). There is evidence of substantial recent genetic yield progress for winter wheat, spring barley, starch potato and sugar beet in the Netherlands (Rijk et al., 2013) and elsewhere in Northwest Europe (Mackay et al., 2011;Peltonen-Sainio et al., 2009). This means the yield gaps estimated by Silva et al. (2017a) were slightly underestimated for some crops and that crop models need to be re-parametrized if they are to provide reasonable benchmarks for farmer field data.
Ya variability during the period 2015-2017 was smallest for starch potato and winter wheat (coefficient of variation, CV, smaller than 15%), intermediate for sugar beet, spring barley and ware potato (CV of 15-20%) and greatest for seed potato and spring onion (CV of 20-25%). The small yield variability of starch potato was expected given that most fields of this crop were located in a single region (northeast) and cultivated in relatively narrow crop rotations (Silva et al., 2017a). Conversely, the small yield variability for winter wheat suggests this is an easier crop to manage across different weather and soil conditions throughout the country. For most crops, Ya were greatest in the polders and in clay rather than sandy soils (Table 2). Water availability during the growing season was also important as summer rainfall was positively associated with the yield of ware potato, onion and wheat. The same was true for spring and summer irrigation in case of ware potato. For some crops, later sowing and earlier harvest were associated with lower yields (Fig. 3). Sowing dates were also identified as an important driver of yield variability in other intensive cropping systems (Rattalino-Edreira et al., 2017). Sowing and harvest dates depend on crop type (Fig. 3), labour and machinery constraints at field/farm level (Reidsma et al., 2015a) and climatic differences between regions across the country. The fact that cereals have a well-defined physiological maturity, while root, tuber or bulb crops do not, allows adjusting the harvest date of the latter to the availability of labour and machinery.
We were not able to analyse the effects of pest, disease and weed management (Lechenet et al., 2017;Andert et al., 2015) or of preceding crop and cropping frequency (Weiser et al., 2018;Silva et al., 2017a;Andert et al., 2016) but expect these to be of key importance when devising strategies to narrow yield gaps in intensive arable cropping systems. Incidence of pests and diseases is the most likely explanation for the relatively large yield gaps of winter wheat while rotational effects, and related soil-borne diseases, play a more prominent role in bulb and tuber crops. As an example, starch potato is often cultivated in the same field once every two to three years and such high cropping frequency is well known to reduce potato yields (Struik and Bonciarelli, 1997;Scholte and s'Jacob, 1990). In summary, there is limited scope to further narrow of yield gaps in Dutch arable cropping systems and interventions to do so should be targeted with caution due to potential environmental externalities. This is especially important due to the small yield responses to nutrients observed for most crops (Table 2; Silva et al., 2017a).

Can water productivity gaps be narrowed?
This study provides water productivity benchmarks for key arable crops in Northwest Europe (Fig. 5). Literature on water productivity for crops such as sugar beet and onion is scant, but we could evaluate our benchmarks for ware potato and winter wheat against those of the Global Yield Gap Atlas (GYGA) and associated references (Rattalino-Edreira et al., 2018;Wang et al., 2018). Our estimates of potential water productivity align well with those from crop model simulations of Yp for ware potato in China (Fig. 7A). For winter wheat, quantile regressions fitted to the 90th percentile of the GYGA water-limited yield data of wheat (10-yr average per weather station) had a slope of 34 kg DM ha −1 mm −1 (Rattalino-Edreira et al., 2018; but the plateau was much greater for the subset of GYGA data for Northwest Europe than for the full GYGA dataset - Fig. 7B). This slope of 34 kg DM ha −1 mm −1 for rainfed wheat is closer to the potential water productivity averaged across fields and years estimated in our study assuming no water limitations during the growing season (ca. 36 kg DM ha −1 mm −1 ; yellow star symbol in Fig. 7B) than the yearspecific slope of the fitted quantile regression (42 kg DM ha −1 mm −1 ; solid line in Fig. 7B). Finally, these potential water productivity values are all much greater than those reported for winter wheat elsewhere (Lollato et al., 2017), which confirms the high potential water productivity for this crop in Northwest Europe.
Our results suggest there is some scope for increasing water productivity of arable crops in the Netherlands (Fig. 5A). Although the total amount of water available during the growing season is enough to avoid major water limitations for most crops, rainfall distribution and availability in key periods of the growing season remain problematic and require adaptation measures at field, farm and regional levels (van Duinen et al., 2015). In the short-term, water productivity can be best achieved by narrowing yield gaps through adoption of precision agriculture practices and provision of supplementary irrigation to avoid water limitations in key periods of the growing season (Silva et al., 2017a). In the long-term, it is important to devise strategies to manage and use the current water surplus in an efficient way. These include improving the current drainage systems in clay soils and improving the capacity to store water, and the efficiency of current irrigation systems, in sandy soils (van Duinen et al., 2015). Although this will all become more important in the context of future climate change, which is expected to further increase rainfall variability in the Netherlands (Reidsma et al., 2015b;Schaap et al., 2013), such interventions require large investments and coordination among stakeholders from farm to regional levels. The ability to narrow water productivity gaps for arable crops in the Netherlands is thus questionable, also given the small scope to further narrow yield gaps.

Opportunities and uncertainties of 'big data'
'Big data' from commercial farms offer new opportunities for agronomic research and to improve farm management (Wolfert et al., 2017). In this study, we coupled a large database of farmer field data to daily weather data and used a summary crop model to derive benchmarks for crop yields and water productivity. Our experiences indicate that at this moment 'big data' are particularly useful to describe current cropping systems in terms of input-output coefficients and to derive benchmarks for farm performance in relation to crop yields (Fig. 3), water productivity (Fig. 5) or N use efficiency (Quemada et al., 2020). However, the relatively low R 2 of the fitted regressions (Table 2) indicates biophysical, variety and management information explained only up to 20-49% of the variation in reported Ya and that such methods are unsuitable for predicting Ya based on this type of information. This points to data quality problems, missing data and/or to the inability of empirical methods to reliably capture the relative contribution of genotype, environmental and management factors to crop yields. The database used in this report contains uncertainties as it builds upon the willingness of farmers to record crop management operations in detail. Data quality issues associated with farmer reported data are well-known (e.g., Fraval et al., 2019;Weiser et al., 2018). In our case, we expect most problems on the recording of fertilizer data (in particular organic sources) and hence, excluded from the analysis all potato and wheat fields without records of organic N applied. Irrigation data may also not be recorded in detail by farmers due to lack of incentives to do so (e.g., low cost of water). We also acknowledge uncertainties in haulm killing dates for potato crops, particularly for seed potato which is killed early in the season (July -August) to ensure the right tuber size and healthy seed. Despite these issues, spatially explicit databases with detailed biophysical (Figs. 4 and 5) and crop management information (Table 2), such as the one used in this study, can be used to describe and derive benchmarks for resource use efficiencies in farmers' fields.
Data availability determines the indicators that can be calculated and their uncertainty. The crop modelling framework used to simulate Yp and the components of the water balance have a number of limitations as documented elsewhere (Allen et al., 1998). However, this pragmatic but widely used approach (e.g., Paredes et al., 2018) was preferred over more mechanistic crop models because the latter are not readily available for all crops studied (e.g., onion) and have not been thoroughly tested against recent field data (Seidel et al., 2018), especially for crops like potato or sugar beet. Further attention should also be paid to the simulation of water-limited yields and the effects of capillary raise, particularly in regions with shallow groundwater tables (Kroes et al., 2018). As such, we think the current approach provides a good starting point to identify the bottlenecks and sources of uncertainty of current benchmarks.

Conclusion
Yield gaps and water productivity are key indicators to monitor progress towards more sustainable cropping systems and large datasets of reported crop yields and management details ('big data') can play a major role in this. This study provides quantitative insights into the magnitude and drivers of these indicators for the main arable crops cultivated in Northwest Europe, building upon individual farmer field data for the Netherlands. Yield gaps were ca. 10% of the potential yield (Yp) for sugar beet, 25-30% of Yp for ware, seed and starch potato and spring barley, and 35-40% of Yp for starch potato, spring onion and winter wheat. In general, further narrowing yield gaps requires finetuning of sowing and harvest dates and coping with the risk of drought in key periods of the growing season, which may not always be possible due to climatic and logistical constraints or environmental legislation. Although not covered in this study, we acknowledge that yield reduction due to pests, diseases, weeds and narrow crop rotations are also key explaining factors to Ya in the intensive cropping systems studied. The omission of these factors may partly explain the relatively low R 2 of the fitted regressions, an issue that merits further research. Results suggest there is some scope to improve water productivity as current values of 13 kg DM ha −1 mm −1 for spring barley, ca. 15 kg DM ha −1 mm −1 for seed potato, spring onion and winter wheat, 23 kg DM ha −1 mm −1 for ware potato and ca. 25 kg DM ha −1 mm −1 for starch potato and sugar beet are about half of their potential. This could be best achieved in the short-term through yield gap closure but in the long-term it is important to devise strategies to store and use surplus water and cope with rainfall variability. Further research should focus on understanding in greater detail the role of climatic extremes, growth-reducing factors and rotational effects on crop yields (Silva et al., 2017a;Schaap et al., 2011) and on balancing current productivity levels with increases in resource use efficiency and decreases in emissions per unit area. Finally, we conclude 'big data' are useful to characterize cropping systems at regional scale and to develop benchmarks for farm performance but not as much to explain yield variability or make predictions in time and space. Fig. 7. Simulated Yp and evapotranspiration estimated in this study vis-á-vis those available in the Global Yield Gap Atlas (www.yieldgap.org) for (A) yearly-specific Yp of potato in China (Wang et al., 2018) and (B) Yw of wheat averaged per weather station over a 10-year period in different world regions (Rattalino-Edreira et al., 2018). Solid lines are as per Fig. 5B and G. The dashed and dotted lines in (B) are fitted quantile regressions to data from NW Europe only and to the full GYGA dataset, respectively.