Hydrological Functioning of Maize Crops in Southwest France Using Eddy Covariance Measurements and a Land Surface Model

Abstract

The primary objective of this study is to evaluate the representation of the energy budget for irrigated maize crops in soil–vegetation–atmosphere transfer (SVAT) models. To this end, a comparison between the original version of the interactions between the soil–biosphere–atmosphere (ISBA) model based on a single-surface energy balance and the new ISBA-multi-energy balance (ISBA-MEB) option was carried out. The second objective is to analyze the intra- and inter-seasonal variability of the crop water budget by implementing ISBA and ISBA-MEB over six irrigated maize seasons between 2008 and 2019 in Lamasquère, southwest France. Seasonal dynamics of the convective fluxes were properly reproduced by both models with R² ranging between 0.66 and 0.80 (RMSE less than 59 W m⁻²) for the sensible heat flux and between 0.77 and 0.88 (RMSE less than 59 W m⁻²) for the latent heat flux. Statistical metrics also showed that over the six crop seasons, for the turbulent fluxes, ISBA-MEB was consistently in better agreement with the in situ measurements with RMSE 8–30% lower than ISBA, particularly when the canopy was heterogeneous. The ability of both models to partition the evapotranspiration (ET) term between soil evaporation and plant transpiration was also acceptable as transpiration predictions compared very well with the available sap flow measurements during the summer of 2015; (ISBA-MEB had slightly better statistics than ISBA with R² of 0.91 and a RMSE value of 0.07 mm h⁻¹). Finally, the results from the analysis of the inter-annual variability of the crop water budget can be summarized as follows: (1) The partitioning of the ET revealed a strong year-to-year variability with transpiration ranging between 40% and 67% of total ET, while soil evaporation was dominant in 2008 and 2010 due to the late and poor canopy development; (2) drainage losses are close to null because of an impervious layer at 60 cm depth; and (3) this very specific condition limited the inter-annual variability of irrigation scheduling as crops can always extract water that is stored in the root zone.

1. Introduction

As the largest consumer of water [1], the agricultural sector is already being threatened by global changes that could increase the frequency of heat-stress [2], drought [3], and water scarcity [4,5]; while food demand rises as well [6]. In France, maize represents only about 17% of cereal production, but uses between 70–80% of irrigation water [7]. This high water demand unfortunately coincides with the period when water demand peaks in the other sectors, while river flows are at their lowest annual levels. Within this context, improving our knowledge of the different water budget terms of irrigated maize including evapotranspiration (ET) is a first step toward the conception of resilient agricultural practices to improve the management of water through better irrigation scheduling and toward future projection of maize water needs. In addition, water is a strategic transverse issue of several sustainable development goals of the United Nations, including goal 2 “zero hunger”, goal 6 “clean water and sanitation”, and goal 13 through the impact of climate change on water resources.

The hydrological functioning of crops, which defines the dynamics of the different terms of the water budget including evapotranspiration and its partition between soil evaporation, plant transpiration, and drainage fluxes can be analyzed using land surface models (LSMs) of various complexities [8,9,10]. The FAO-56 double-coefficient approach [11] solves a simple soil water budget and estimates ET by adjusting the climatic evaporative demand by an empirical coefficient. The latter called “the cultural coefficient” encompasses all the processes related to crop growth that are not explicitly represented by the model. This approach has been extensively assessed by comparison with in situ measurements [12,13], including over maize crop [14] and used for water management and irrigation scheduling purposes [15,16]. While an increased referenced-evaporative demand has already been projected at various horizons under different future scenarios [17], projecting the cultural coefficients under future climate conditions is still an open question.

Beside these conceptual approaches, other appealing tools to compute ET fluxes are the surface energy budget models driven by remotely sensed land surface temperature (LST) [18,19,20]. With regards to the above-mentioned water balance approaches, irrigation water inputs, which are very uncertain at the field scale, are not required, as the LST is a proxy of the crop hydric status when water is limited [21,22,23]. Indeed, energy budget models provide an instantaneous estimate of ET at the time of the satellite overpass meaning that (1) an extrapolation between two satellite acquisitions is needed [24] and (2) it is not possible to project ET under future climate change scenarios. In addition, there are no available LST observations combining high spatial resolution (<100 m) and high revisit time (<5 days) with the sensors currently in orbit. Finally, the physically based soil–vegetation–atmosphere transfer (SVAT) models that were initially developed to provide atmospheric models with the lower boundary conditions for weather forecasting or climate projection applications jointly solve the water and the energy budgets of the land surfaces [25,26,27,28]. While rich literature dealing with the evaluation of SVAT models for natural ecosystems exists [29,30,31], there is still a need for assessment studies of such modelling tools over agro-systems characterized by highly transient moisture conditions and a rapid switch of energy between sensible and latent heat fluxes associated with irrigation practices [32,33]. This is probably due to the limited number of agro-system sites within the networks of eddy covariance towers, e.g., Ameriflux [34], Ozflux [35], and FLUXNET [36]. Nevertheless, improving the representation of agro-systems in land surface models is becoming a key issue to improve future climate projection [37].

The land surface processes within SVAT models encompasses a large variety of parameterizations for representing the way vegetation controls water and energy exchanges within the soil–vegetation–atmosphere continuum [26,38,39]. However, within these models, the structural representation of the surface energy budget is a key issue [40,41]. This structural representation involves how the canopy and the soil sources are being coupled and how their exchanges with the atmosphere are being treated. A large part of SVAT models solve one composite energy budget for soil and vegetation such as the interactions between the soil–biosphere–atmosphere (ISBA) [27]; and the Tiled ECMWF Scheme for Surface Exchanges over Land (TESSEL) [42]; meaning that the same temperature is attributed to both components. This “single” source approach faces several limitations for heterogeneous canopies, with part of the bare soil directly exposed to incoming radiations [41,43] when contrasted hydric conditions exist between a dry surface layer and a well-watered root-zone layer, which is a very likely situation for irrigated crops. In this case, a strong difference of temperature is expected between the dry, the hot soil surface, and the transpiring vegetation. This impacts the partitioning of the net radiation between latent and sensible heat fluxes [18] and the partitioning of ET between soil evaporation and plant transpiration. Within this context, several models solving two separate energy budget for soil and vegetation were developed such as the new multi-energy balance version of ISBA named ISBA-MEB [44] or the simple biosphere model [45]. For irrigated maize crops, the choice between a single- and a double-source energy budget approach is still an open question.

Within this context, the objective of this study is two-fold:

• Assessing the ability of an SVAT model to reproduce the different terms of the energy and water budget over six maize seasons. To this objective, the original version of the SVAT model named interactions between the soil–biosphere–atmosphere (ISBA) [27] that solves a single energy budget is compared with the new multi-energy balance version (ISBA-MEB). ISBA-MEB provides a good environment to compare single- and dual-source representations of irrigated maize as all other processes are parameterized in the same way in both versions of the model.
• Investigating the year-to-year variability of the different terms of an irrigated maize with a special emphasis on the water losses for the plant (drainage and soil evaporation).

2. Materials and Methods

2.1. Experimental Site

Lamasquère (FR-Lam) (43°50′05″ N, 1°24′19″ E) is a well-instrumented plot owned by the Purpan Engineering School (Toulouse, France) and managed as a milk production farm. This 24 ha site (Figure 1) is located in the southwest of France (Haute Garonne department), and it is a class 1 site within the Integrated Carbon Observation System (ICOS) network (https://www.icos-cp.eu/; (accessed on 23 May 2021)) and regional spatial observatory (OSR) project [46]. Positioned along the Touch river, it is an intensively cultivated site dominated by western and southeastern winds, with an average wind speed of 1.8 m s⁻¹ over the site [47].

Over the current study period, FR-Lam experienced a mean annual precipitation of 642 mm and a mean air temperature of 12.8 °C. The soil texture consists of loam (33.7%), clay (54.3%), and sand (12%) [48], and it is vertically heterogeneous over the depths at which we had observations. There is a strong crop rotation on this field, and only cropping years with maize were selected (2008, 2010, 2012, 2014, 2015, and 2019). This summer crop is usually sown between April and early May. Irrigation is customarily applied to substitute the depleting soil water from early June, and it takes about 4–5 days to cover the full plot. The pivot irrigation method was applied in 2008 and 2010, while for the other years, the farmer adopted the sprinkler technique. At the commencement of senescence, harvest often occurs between August and September of the same planting year for silage. A crop season in this study is taken from 1st April to 30th September.

Table 1. Maize crop cycles, total amount of rainfall and irrigation, sowing and harvest dates, maximum Leaf Area Index (LAI).

Crop Season	Rain (mm)	Irrig (mm)	Sowing Date	Harvest Date	Max LAI (m2 m−2)
01/04/2008–30/09/2008	361.6	27.3	20/05/2008	12/09/2008	3.89
01/04/2010–30/09/2010	277.6	106.1	21/04/2010	20/09/2010	4.01
01/04/2012–30/09/2012	292.5	148.7	27/04/2012	27/08/2012	5.89
01/04/2014–30/09/2014	256.2	167.7	14/05/2014	24/09/2014	5.21
01/04/2015–30/09/2015	261.5	149.5	05/05/2015	08/09/2015	6.64
01/04/2019–30/09/2019	267.4	152.4	22/04/2019	12/09/2019	4.8

highlights the field operations and the maximum leaf area index (LAI) recorded for the selected crop seasons.

2.2. Data Set Description

2.2.1. Eddy Covariance and Meteorological Measurements

At the center of FR-Lam (except for 2019), an Eddy Covariance (EC) system is installed at a height of 3.65 m. This system continuously monitors turbulent fluxes combining synchronized 3D wind components from a 3D sonic anemometer (Campbell, CSAT 3), and the turbulent component of the scalar of interest: air temperature using the same instrument and water vapor concentration from the open-path infrared gas analyzer (LiCor Li-7500, IRGA). These measurements were carried out within the footprint of the site oriented in the direction of the prevailing wind (west and southeast). Half-hourly sensible heat (H, W m⁻²) and latent heat fluxes (LE, W m⁻²) were estimated with the EDIRE software in accordance with the CarboEurope-IP protocol [47]. The flux measurements for 2019 were however obtained with the standard ICOS set-up using both the closed path IRGA (Li-7200) from LiCor and a 3D-sonic anemometer (HS50) from Gill Ltd. installed on a 5 m mast due to a taller maize species. Missing fluxes were gap-filled using marginal distribution sampling according to [49] and [50]. These half-hourly flux data were corrected in accordance with the statistical and objective criteria described in [51].

Meteorological sensors measured the net radiation (CNR1, Kipp and Zonen), wind speed and direction (Wind-vane/prop Young), air temperature, relative humidity (HMP35, Vaisala), and precipitation at a half-hourly time step. Erroneous and missing data were reconstructed in accordance with the approach in [51].

Volumetric water content values were also collected at the following depths (0, 5, 10, 30, 50, and 100 cm) by a TDR CS616 from 2005 to 2011 and by a thetaprobe ML2X from 2012 to 2019 with 4 replicates located around the mast (see Figure 1). These sensors were calibrated using the gravimetric method. In addition, an electric piezometer (NIVO-05 model) with a pressure transducer device was taking measurements of both the pore pressure and the ground water level. Manual measurement of the water table is performed fortnightly to check for drifts, and it is recalibrated when necessary. For further information on the instrumentation of this site, see [47].

2.2.2. Ground Heat Flux

The ground heat flux (G) was taken as the average of measurements from four HFPs (see B, C, D, and E in Figure 1) in order to guarantee representativeness of the source area. Furthermore, using the calorimetry approach described in [52], G was corrected with the energy stored between the soil surface and HFP.

Due to instrumental failures and field operations, there were several missing data in 2014, 2015, and 2019. The gap-filling of G was carried out as follows: (i) When the LAI is less than 0.5 m² m⁻², G was set as a constant fraction of the net Radiation (Rn) following [53] (see Equation (1a)). (ii) When LAI is greater or equal to 0.5 m² m⁻², the empirical method developed by [54] which accommodates the shielding effect of the canopies on the soil (Equation (1b)) was used.

$$G=a_1{R}_n \mathrm{if} \mathrm{LAI}<0.5 {\mathrm{m}}^2 {\mathrm{m}}^{-2}$$

(1a)

$$G=[a_3+a_3{e}^{-a_{4}LAI}] R_n \mathrm{if} \mathrm{LAI}\ge 0.5 {\mathrm{m}}^2 {\mathrm{m}}^{-2}$$

(1b)

The calibration of these two empirical expressions for the crop years that have sufficient G measurements (2008, 2010, and 2012) was done by splitting the data set into a 70% training set and a 30% testing set using a Jacobian-based optimization method to minimize the root mean square error (RMSE).

For the low LAI category, the optimized coefficient is a₁ = 0.29. This is close to the value used in the literature by different models [55]. This value resulted in an R² of 0.80, 0.71, 0.78, and 0.76 for the years 2008, 2010, 2012, and the combined years, respectively. In addition, low RMSE values of 15, 18, 15, and 16 W m⁻² were obtained, respectively.

Likewise, in the second category of high LAI, the optimal coefficients (a₂ = 0.049; a₃ = 0.135; a₄ = 0.249) resulted in R² of 0.89, 0.76, 0.88, and 0.83; RMSE of 9,14, 9, and 11 W m⁻² for 2008, 2010, 2012, and the combined years, respectively.

2.2.3. Sap Flow Measurements

Sap flow measurements are used in this study to assess the performance of the models in the partitioning of ET into soil evaporation—the water lost and plant transpiration- the water used by the crop. Twenty maize stands (all located within the footprint of the EC) were equipped with micro sapflow sensors (Dynamax Inc.’s SGB19 and SGB25) from 2 July 2015 to 23 August 2015 when the plants were close to physiological maturity (min LAI = 4.0 m² m⁻², min crop height = 1.59 m, mean diameter = 22.4 ± 3.4 mm). The selected maize stands were assumed to be representative of the field. Measurements were performed using the Heat Balance Method (HBM) that was developed by [56] and [57] based on a thermal flow meter approach. To extrapolate sap flow measurements over the entire plot for comparison with the EC measurements, transpiration values were obtained by multiplying the averaged sap flow for the 17 retained sensors by the plant density (8.6 plants m⁻²). Readers are referred to Figure S1 in the Supplementary Materials for the time series of the sap flow data.

2.2.4. Vegetation Characteristics and Irrigation Water Amounts

During each growing season, the vegetation LAI was monitored about five times with higher frequency during the rapid growth phase. The sampling process involved collecting vegetation on 10–20 subplots, and then the LAI is measured using a planimeter (Li-3100, LiCor, Lincoln, NE, USA) [58]. More details on the protocol can be found in [58]. The LAI, vegetation height (hveg), and fraction cover (fc) measurements were linearly interpolated with a 10 days’ time step. Irrigation quantity and timing were also gathered through regular surveys.

2.2.5. Water Storage Calculation

Equation (2) highlights the principal components of a hydrological budget. D and R represent the drainage and surface run-off, respectively [59]. Using Equation (3), the soil water storage (S) was estimated as an integration between 0 and 50 cm (30 cm in 2008 and 2010 due to missing data), since most of the water taken up by maize is between these depths [60]. θ is the soil water content in m³ m⁻³ at layer i, and Z is the layer thickness between layer i and the preceding layer in meters. In 2019, and during the summer of 2014, S is computed using simulated soil moisture (ISBA-MEB) as a result of large gaps in the field measurements. In this study, the fluctuation in the soil water storage (ΔS) was computed as the difference between the storage at the beginning of the crop season and the storage value a few days before harvest due to unavailable soil water data (removal of buried sensors for harvest).

P + Irrig + ΔS = ET + D + R

(2)

$$\mathrm{S} ={{\displaystyle \sum }}_{i=1}^5{\theta }_i{Z}_i^{\prime }$$

(3)

Given the brevity of the irrigation events, the concise amount, and the flatness of the site, R is considered negligible in this study. Hence, D is the residual of the water budget.

2.2.6. Definition of Metrics

The performance of both ISBA and ISBA-MEB were assessed in terms of the RMSE, the coefficient of determination (R²), and the mean absolute error (MAE) as shown in Equation (4a–c):

(i)

RMSE: Represents the root mean square deviation between the measured and the simulated variable.

(ii)

R²: Expresses the proportion of variance in the simulated variable that is predicted from the observed.

(iii)

MAE: The average magnitude of errors between the observation and the simulated datasets.

$$RMSE=\frac{\sqrt{{{\displaystyle \sum }}_{i=1}^N{\left(si{m}_i-ob{s}_i\right)}^2}}{N}$$

(4a)

$$R^2=\frac{n{{\displaystyle \sum }}_{i=1}^N\left(si{m}_{i}ob{s}_i\right)-({{\displaystyle \sum }}_i^{N}si{m}_i)({{\displaystyle \sum }}_i^{N}ob{s}_i) }{\sqrt{\left[N {{\displaystyle \sum }}_{i=1}^{N}si{m}_i^2-{({{\displaystyle \sum }}_{i=1}^{N}si{m}_i)}^2\right]} \sqrt{\left[N {{\displaystyle \sum }}_{i=1}^{N}ob{s}_i^2-{({{\displaystyle \sum }}_{i=1}^{N}ob{s}_i)}^2\right]}}$$

(4b)

$$MAE=\frac{{{\displaystyle \sum }}_{i=1}^N\left|si{m}_i-ob{s}_i\right|}{N}$$

(4c)

where sim and obs are the simulated and observed values, respectively, N is the number of the observation, and i represents the index of a given variable.

3. Model Description and Implementation

3.1. Model Description

The ISBA model, an embedded component on the SURFEX (in French “Surface Externalisée”) platform [61] is a composite soil–vegetation land surface scheme which simulates surface energy exchanges within the soil–vegetation–atmosphere continuum. In this study, the ISBA-A-gs scheme, which is the ISBA model coupled with the stomatal conductance (A-gs) model of [62,63], was used. ISBA-A-gs considers the functional relationship between atmospheric carbon concentration (CO₂) and stomatal aperture at a leaf level. It also employs the physiological parameterization of the stomatal conductance (g_S) as proposed in [64] and [65] to describe plants’ photosynthesis. Likewise, it accounts for the response of the plants to soil–water stress by employing a normalized soil moisture factor to g_S. In this study, the non-interactive vegetation option is chosen—meaning that the vegetation characteristics (LAI, fc and hveg) are prescribed from in situ measurements. The solar radiation transfer scheme uses a multi-layer configuration that considers sunlit and shaded leaves. The soil’s transfers and heat exchanges are simulated using the recent multi-layer soil diffusion scheme developed by [66,67] which was shown to perform better on crops than the initial force-restore approach [68] by [69]. This is because the force restore method does not account for time-variability of the plants root depth, and it represents the soil using a limited number of layers [69]. The soil grid is composed of 8 discrete soil layers (0, 5 cm, 10 cm, 30 cm, 50 cm, 1 m, 1.5 m, and 2 m) which correspond to the depths with available in situ measurements (0 cm to 1 m). For an extensive description of the model, the governing equations, and the algorithm of the A-gs model, the reader is referred to [62,64,65]. For ease, ISBA-A-gs would be referenced as ISBA henceforth.

The multi-energy balance (MEB) option referred to as ISBA-MEB in this study consists of a multi-source approach developed by [44] and [70] that solves the energy budget for the soil and vegetation separately following [71]. This approach introduces an additional prognostic variable (vegetation temperature), while other forcing variables and processes remain the same as in ISBA.

3.2. Model Implementation

Forcing Parameter and Data

The main input parameters and variables of the models are shown in

Table 2. Values and ranges of the input variables for ISBA and ISBA-MEB using the Ag-s photosynthesis option.

Symbol	Description	Range	Unit	References
fc	Vegetation fraction cover	0–1.0	-	Estimated
LAI	Leaf area index	0–6.6	m2 m−2	Measured
hveg	Vegetation height	0–3.3	m	measured
Emissoil	Soil emissivity	0.90	-	-
Emisveg	Vegetation emissivity	0.97	-	-
albsoil	Soil albedo	0.12–0.16	-	Derived from measurement
albveg	Vegetation albedo	0.17–0.28	-	Derived from measurement
Cl	Sand fraction	0.10–0.12	-	measured
Sd	Clay fraction	0.48–0.53	-	measured
gm	Mesophyll conductance	0.0020	m s−1	[63]; Table 1
Γ	Coefficient for the calculation of the surface stomatal resistance	0.02		[63]; Table 1
Wrmax	Coefficient for maximum interception water storage capacity	0.2		[27]; Equation (24)
Cv	Thermal coefficient for the vegetation canopy	0.00005	K m2 J−1	[27]; Equation (8)
gc	Cuticular conductance	0.0004	m s−1	[73]
Θc	Critical normalized soil water content for stress parameterization	0.334	-	[63]; Equation (9)
Dmax	Maximum air saturation deficit	0.065	kg kg−1	[63]; Equation (A3)

. ISBA and ISBA-MEB permit the use of 4–19 patches representing different vegetation types, while bare soil, rocks, and snow are represented by patch 1–3, respectively [61]. Each vegetation patch is associated with a set of parameters that is particularly related to the photosynthesis module provided by the ECOCLIMAP II database [72]. Patch 8, which corresponds to C4 crops, has been used in this study. Similarly, both models were forced with gap-filled half-hourly meteorological variables of air temperature (TA in K), specific humidity (QA in kg kg⁻¹), air pressure (PS in Pa), longwave radiation (LWin W m⁻²), wind speed (WS in m s⁻¹), rainfall (RAIN in kg m⁻²), and shortwave radiation (SWin W m⁻²).

Total surface albedo ($al{b}_{total}$) was estimated as an instantaneous ratio of up-welling radiation to down-welling radiation measured by the CNR4 net radiometer between 10 h 30 and 15 h to avoid low solar zenith angles as indicated in [74]. Afterwards, albedo values for soil ($al{b}_{soil}$) and vegetation ($al{b}_{veg}$) were obtained using Equation (5). The same distinction was applied to emissivity. For each crop season, this study kept alb_soil and alb_veg constant independent of soil moisture and vegetation development, respectively.

$$al{b}_{total}=[f_{c}al{b}_{veg}+\left(1-f_c\right)al{b}_{soil}]$$

(5)

The pedotransfer function often used in ISBA [75] resulted in a large under-estimation of the transpiration term as already pointed out by [69]. Hence, the soil hydraulic properties were computed using the pedotransfer function based on the measurements of the sand and clay content at the six depths [76].

Table 3. Soil texture and soil hydraulic properties for the sap flow period.

Layer (cm)	Clay (%)	Sand (%)	CONDSAT	MPOTSAT	BCOEF	WP (m3 m−3)	FC (m3 m−3)	SAT (m3 m−3)
0	51	12	0.0000001081	−0.5551	10.488	0.1052	0.2520	0.2913
5	52	11	0.0000001067	−0.5665	10.625	0.1062	0.2550	0.2924
10	53	10	0.0000001057	−0.5781	10.762	0.1152	0.2556	0.2933
30	50	11	0.0000001088	−0.5665	10.351	0.1045	0.2420	0.2933
50	48	11	0.0000001120	−0.5665	10.077	0.1033	0.2289	0.4222

presents a summary of the hydraulic properties used during the sap flow period. The initial soil water values at the start of the season are forced using in situ measurements at each layer.

Similarly, a modification of the model was done by modeling a shallow saturated zone [77], which accommodated the ground water reservoir at the site. Initial runs of both models revealed saturated zones at shallow depths, indicating a water table (WT) that is shallower than the lower boundary condition used in ISBA. The high water retention capacity of the soil in FR-Lam, coupled with the nearby presence of the Touch river, encourages the existence of such a shallow WT which counters the drainage lower-boundary condition that ISBA-DIF employs. For this reason, the measured WT was used as a lower-boundary condition. This approach provided a more realistic estimation of both water and surface energy fluxes due to the optimal capturing of capillary processes. See [77] for further information on this method.

The maximum value for the root depth was set to 1.5 m following [78], and 87% of the root was apportioned between 0 cm and 50 cm following [79], while the other fraction was allocated to the deeper layers.

4. Results and Discussion

4.1. Experimental Data Analysis

4.1.1. Meteorological Conditions and Vegetation Characteristics

Figure 2 provides a general overview of the climate and the crop growth conditions along the growing seasons. There is a strong seasonality of air temperature. In addition, while the inter-annual variability of rainfall (+ irrigation) is small (from 388 mm in 2010 to 441 mm in 2012), the seasonal variability of the irrigation is significant. The highest irrigation event occurs around the peak of the crops’ growth, between July and September concomitantly with the highest temperature. The farmer applied no irrigation in June for 2010 and 2012 thanks to abundant rainfall during the spring months in contrast to the other growing seasons. The year 2019 was characterized by a drought period particularly in July and exhibited the highest cumulative water inputs during July and August of all the growing seasons. Likewise, a strong inter-annual variability of the LAI cycles is observed with a rapid growth in 2015 thanks to favorable temperature conditions in April (13.2 °C) and in May (16.2 °C), while a late growth is evident in 2010 mainly because of the sub-normal air temperature (see Figure 2b). In 2008, the crop may have suffered from either a very late sowing, or the cultivation of a different species. The peak LAI values also differ strongly from one year to another with exceptional growing conditions in 2015 leading to a peak LAI of 6.6 m² m⁻² thanks to high irrigation inputs applied in July by contrast with 2008 and 2010 characterized by LAI below 4.1 m² m⁻². Finally, a significant shift up to 1 month is also observed for the harvest date

4.1.2. Energy Balance Closure

Numerous studies on surface energy balance have clearly hinted a non-closure due to several reasons [80], including in our study site [81]. Before comparing in situ measurements with simulated values, a quality check of the EC flux was carried out using the energy balance closure (EBC) as an indicator. This proxy is the slope obtained from the statistical linear regression of (LE+H) against (Rn-G). Figure 3a shows the monthly closure (box plot) and residual energy (black circles) from April to September. The spring months had the lowest mean closure (median of 0.61 in April), while a better closure is observed in summer (0.73 in August) when vegetation is fully developed. This seasonality effect has been described in [81]. Mean monthly residual energy varied between 37 and 51 W m⁻² with the smallest and largest residual energy found in 2015 (22 W m⁻²) and 2008 (70 W m⁻²), respectively.

These statistical results are within the reported values associated with agricultural sites in this region [82,83], and this shows good data quality from the EC system. Before making comparisons with the estimations from the models, the Bowen ratio correction method was applied to the raw fluxes [84]. This method assumes that the Bowen ratio (β) is correctly measured by the EC system and partitions the deficit energy into both H and LE. This method has been successfully applied in numerous studies [85].

Figure 3b shows the closure and the residual energy after correction. Table S1 in the Supplementary Materials summarizes the scores of the raw fluxes (raw) in comparison with the corrected fluxes (corr) at a half-hourly time step. As expected, the EBC of the corrected fluxes is higher than 0.90, and the residual energy was not greater than 14 W m⁻². In subsequent sections, these corrected fluxes would be referred to as the in situ measurements.

4.2. Assessment of the ISBA and ISBA-MEB Models

4.2.1. Energy Budget

The 2-days average time series (between 12 h and 15 h 30) comparing the measured and simulated fluxes of Rn, LE, H, and G are displayed in Figure 4, Figure 5, Figure 6 and Figure 7, respectively for the 6 crop seasons. The statistical metrics at half-hourly time steps are reported in

Table 4. Statistical metrics comparing the energy budget components estimated by ISBA and ISBA-MEB against in situ measurements (units of RMSE and MAE in W m⁻²) at half-hourly time step. In bold fonts are cases when ISBA-MEB gave better scores (R²).

Maize Year	Model	Indicator	Rn	LE	H	G
2008	ISBA	R2	0.97	0.81	0.67	0.75
		RMSE	34.0	56.0	42.9	37.1
		MAE	25.0	33.9	27.5	26.6
	ISBA-MEB	R2	0.98	0.83	0.73	0.81
		RMSE	33.5	50.7	30.2	32.1
		MAE	27.5	31.7	25.9	27.1
2010	ISBA	R2	0.98	0.81	0.73	0.81
		RMSE	33.4	46.9	45.1	28.8
		MAE	22.1	31.8	29.9	21.5
	ISBA-MEB	R2	0.98	0.84	0.77	0.89
		RMSE	32.2	42.2	31.2	21.7
		MAE	22.6	29.7	27.1	26.9
2012	ISBA	R2	0.96	0.77	0.66	0.74
		RMSE	42.7	58.9	58.5	26.7
		MAE	27.5	36.1	39.3	34.7
	ISBA-MEB	R2	0.96	0.80	0.68	0.86
		RMSE	41.7	53.2	43.3	19.8
		MAE	28.2	34.0	34.3	43.2
2014	ISBA	R2	0.96	0.79	0.78	0.73
		RMSE	46.6	50.9	47.4	38.1
		MAE	28.3	31.5	34.5	29.7
	ISBA-MEB	R2	0.96	0.85	0.80	0.79
		RMSE	46.7	42.7	37.6	41.1
		MAE	29.4	28.1	30.6	32.9
2015	ISBA	R2	0.97	0.87	0.70	0.76
		RMSE	40.6	43.1	51.8	41.9
		MAE	29.5	34.1	45.8	39.7
	ISBA-MEB	R2	0.97	0.88	0.72	0.79
		RMSE	39.7	39.6	41.8	44.8
		MAE	29.9	33.0	40.9	48.5
2019	ISBA	R2	0.96	0.85	0.76	0.76
		RMSE	44.2	44.3	47.6	37.4
		MAE	30.1	26.7	33.7	36.7
	ISBA-MEB	R2	0.96	0.88	0.78	0.82
		RMSE	43.0	38.9	34.2	41.5
		MAE	31.7	26.3	28.1	48.3

.

• Net radiation

Both models presented a good agreement in magnitude with the observations, and the seasonal dynamics were properly reproduced (Figure 4) with R² ranging between 0.96 and 0.98, while the MAE remained lower than 32 W m⁻² (Table 4)_. This was expected as albedo is the main governing parameter, and it has been calibrated based on the observations of the incoming and outgoing shortwave radiation in this study. In addition, metrics of both models are close; although, ISBA-MEB provided a slightly lower RMSE on average. Considering that the shortwave component of Rn is computed similarly in both models, this observed difference arises from the outgoing longwave component seen in the difference between the composite (surface) temperature predicted by ISBA and the explicit (soil and vegetation) temperatures predicted by ISBA-MEB (not shown).

2.

Latent heat and sensible heat flux

In the early crop stages (April and May), LE fluxes show weak values in response to the magnitude of Rn as shown in Figure 4 and Figure 5. Afterwards, and up until harvest, a minimum of 70 W m⁻² is maintained due to the presence of the transpiring plants and highly humid soil. The seasonal dynamics and the magnitude of the in situ measurements are well captured, and the differences between both models are small on average over the season. However, the ISBA-MEB model performed better with R² values between 0.80 and 0.88, while for ISBA, R² was between 0.77 and 0.87 as shown in Table 4.

Interestingly, ISBA was unable to capture the transition period from low crop to developed crop when the heterogeneity of the field was high; this resulted in a large underestimation of LE, with an average MAE of 69 W m⁻² during these time frame in comparison to ISBA-MEB (48 W m⁻²). During these transition periods, average RMSE in ISBA (49 W m⁻²) was 28% higher than that of ISBA-MEB (38 W m⁻²), confirming the need for an explicit representation of the soil and the vegetation. Afterwards, when the field becomes homogenous with the vegetation fully shielding the soil, both models were in better agreement apart from some small discrepancies found during the senescence stage. This is consistent with the findings of [70]; for 87% of the studied forested FLUXNET sites, the difference between the RMSE of ISBA and ISBA-MEB is more significant for sites with low LAI.

Similar conclusions can be drawn from the sensible heat time series (Figure 6), but the difference between both models during the aforementioned transition period is higher than for LE. ISBA predictions of H are characterized by an R² ranging from 0.66 to 0.78, while R² is between 0.68 and 0.80 for ISBA-MEB. Between April and May, when the site is almost bare, both models reproduced H considerably well, likewise, when the vegetation is fully developed. By contrast, during the growth period, a notable discrepancy is observed characterized by an over-estimation of H by ISBA. During these specific periods, ISBA over-predicted by 43, 38, 57, 56, 59, and 41 W m⁻² for 2008, 2010, 2012, 2014, 2015, and 2019, respectively; while the estimations of ISBA-MEB were closer to the in situ measurements with smaller bias (7 and 39 W m⁻²). A similar observation was reported over an open olive orchard by [32], who reported a significant improvement with ISBA-MEB in the prediction of convective fluxes over open canopies characterized by a mix of bare soils directly exposed to incoming radiation and of vegetation shielding the underlying soil. ISBA-MEB being a dual-source model couples both heat sources, thus, better suited during this transition phase.

3.

Ground heat flux

Figure 7 shows that the simulated values of the models are in good agreement with the measurements. The statistics Table 4, taken at a half-hourly resolution, reveals that the accuracy of ISBA is lower than that of ISBA-MEB as R² on average increased from 0.76 to 0.83, while the mean RMSE is reduced from 35 to 33 W m⁻².

4.2.2. Soil Moisture Time Series

Soil moisture is a key variable that integrates all water balance components. For illustration purposes, Figure 8 shows the daily time series of the observed and the predicted soil water content during the sap flow period, and

Table 5. Statistical metrics of the comparison between the observed soil moisture and the estimates of ISBA and ISBA-MEB, over the six maize years. The units of RMSE and MAE are in mm and in bold fonts are cases when ISBA-MEB had better metrics (R²).

Year	Indicator	Model	0 cm	5 cm	10 cm	30 cm	50 cm	100 cm
2008	R2	ISBA	-	0.77	0.80	0.71	-	0.30
		ISBA-MEB	-	0.77	0.79	0.75	-	0.30
	RMSE	ISBA	-	0.024	0.021	0.019	-	0.029
		ISBA-MEB	-	0.026	0.023	0.019	-	0.032
	MAE	ISBA	-	0.026	0.032	0.053	-	0.060
		ISBA-MEB	-	0.029	0.037	0.056	-	0.062
2010	R2	ISBA	-	0.60	0.71	0.46	-	0.40
		ISBA-MEB	-	0.59	0.76	0.51	-	0.39
	RMSE	ISBA	-	0.033	0.024	0.022	-	0.022
		ISBA-MEB	-	0.035	0.024	0.023	-	0.025
	MAE	ISBA	-	0.027	0.019	0.037	-	0.035
		ISBA-MEB	-	0.029	0.020	0.038	-	0.037
2012	R2	ISBA	0.40	0.59	0.40	0.30	0.91	0.32
		ISBA-MEB	0.40	0.63	0.49	0.30	0.89	0.31
	RMSE	ISBA	0.060	0.032	0.032	0.025	0.019	0.027
		ISBA-MEB	0.062	0.034	0.034	0.025	0.023	0.033
	MAE	ISBA	0.059	0.031	0.039	0.025	0.065	0.045
		ISBA-MEB	0.064	0.033	0.04	0.025	0.079	0.052
2014	R2	ISBA	0.42	0.42	0.68	0.33	0.73	0.38
		ISBA-MEB	0.46	0.40	0.76	0.36	0.77	0.47
	RMSE	ISBA	0.046	0.018	0.024	0.015	0.011	0.008
		ISBA-MEB	0.049	0.019	0.023	0.015	0.012	0.008
	MAE	ISBA	0.037	0.027	0.024	0.022	0.027	0.016
		ISBA-MEB	0.039	0.026	0.025	0.022	0.032	0.019
2015	R2	ISBA	0.65	0.57	0.30	0.30	0.48	0.69
		ISBA-MEB	0.71	0.64	0.30	0.26	0.51	0.69
	RMSE	ISBA	0.031	0.021	0.031	0.017	0.013	0.007
		ISBA-MEB	0.023	0.022	0.032	0.029	0.016	0.010
	MAE	ISBA	0.027	0.017	0.045	0.033	0.011	0.013
		ISBA-MEB	0.036	0.018	0.043	0.030	0.015	0.017
2019	R2	ISBA	0.31	0.46	0.53	0.41	0.62	0.58
		ISBA-MEB	0.30	0.45	0.54	0.41	0.67	0.57
	RMSE	ISBA	0.043	0.036	0.014	0.023	0.023	0.012
		ISBA-MEB	0.045	0.038	0.017	0.026	0.043	0.016
	MAE	ISBA	0.046	0.033	0.026	0.023	0.034	0.031
		ISBA-MEB	0.051	0.036	0.024	0.025	0.036	0.036

presents the statistical metrics for the whole study period. Overall, Table 5 showed that: (1) the ability of both models to predict soil moisture is very variable from one year to another; (2) better prediction are obtained with ISBA-MEB; and (3) predictions of both models are acceptable with RMSE remaining below 0.04 m³ m⁻³ without giving heavy considerations to the 0 cm layer because sensors can easily be moved and exposed to the open air.

Figure 8 shows that for all the layers, the dynamics of the simulated soil water responded to every rainfall event. Similarly, a good agreement is observed between the observed and the simulated with a coefficient of determination of 0.81 and 0.77 at the 0 cm and 5 cm depth and even higher at the 50 cm depth (0.84). However, at 10 and 30 cm, a correlation lower than 0.40 was found which might have been affected by the presence of a large clay content just above the 10 cm depth. This limited water infiltration, but the presence of cracks permitted soil–water redistribution. The prediction quality of both models is similar here due to the aforementioned reasons in Section 4.2.2.

4.2.3. Transpiration

Figure 9 presents the daily time series of the transpiration (Figure 9a), the average daily cycles by 14-days period (Figure 9b), and the cumulative values over this period (Figure 9c). When compared with the total ET (not shown), daily ET values lie between 1.6 and 9.5 mm of water while the sap flow values ranged between 0.13 and 7.36 mm. This implies that on average, over this period, transpiration represented 88%, 84%, and 87% of total ET for the in situ data, ISBA, and ISBA-MEB estimates, respectively. These values are close to those obtained in the literatures for a mature irrigated maize with full canopy coverage [86]. These estimates also compare well with the study of [87], which showed that during the vegetation peak (LAI = 4 m² m⁻²), transpiration reached 75% of the total ET for maize cropped in northern China; considering that LAI is above 6 m² m⁻² in our study. As soil moisture strongly influences sap flow, cold fronts indicated by low incoming radiation and cool air temperature (both not shown) resulted in a drop in sap flow in 30 July, 9 August, and 14 August. Likewise, a gentle decline in transpiration is observed in early August, which is influenced by biotic factors such as senescence.

Figure 9b highlights the ability of both models to effectively reproduce the daily dynamics of the sap flow, and its magnitude with no strong time-shift with the observations apart from an earlier drop of transpiration from midday particularly prominent during the first half of August. Both models overestimated transpiration during the first period, and under-estimated afterwards as seen in Figure 9c. Overall, both models underestimated the sap flow measurement by ~24 mm with a slight advantage to ISBA-MEB.

Table 6. Statistical scores on an hourly-scale between the sap flow measurements and transpiration estimated by ISBA and ISBA-MEB from 02/07/2015–27/08/2015.

Model	R2	RMSE (mm h−1)	MAE (mm h−1)
ISBA	0.90	0.074	0.052
ISBA-MEB	0.91	0.071	0.048

shows the statistical metrics of the comparison of the sap flow measurements with the estimations of the models, and it reveals that the predictions of both models are quite similar. This is probably because the sap flow campaign was carried out over a mature and homogeneous maize field (LAI > 4 m² m⁻²). Nevertheless, a slightly better estimation was provided by ISBA-MEB with R² of 0.91 and a lower RMSE of 0.071 mm h⁻¹; whereas, ISBA has R² of 0.90 and RMSE of 0.074 mm h⁻¹. In addition, a bias of 0.052 mm h⁻¹ and 0.048 mm h⁻¹ was obtained for ISBA and ISBA-MEB, respectively. Low transpiration values are also well reproduced by the models, particularly during the first half of August, which is characterized by two periods with a complete stop of transpiration.

4.3. Inter-Annual Variability of Maize Water Budgets

This section assesses the inter-annual variability of the crop water budget, including the partitioning of the ET. Given that ISBA-MEB provided better metrics than ISBA, only the former results are shown in this section.

4.3.1. Evapotranspiration Partitioning

Figure 10 highlights the partitioning of ET (of ISBA-MEB) into its individual components of plant transpiration (E_tr), soil evaporation (E_s), and vegetation evaporation (E_i) at seasonal and monthly time scales.

Figure 10a indicates a strong year-to-year variability of ET partitioning. For example, E_s is the dominating component in 2008 accounting for over 52% of total ET, while E_tr hardly reached 40%. In 2010, a very close contribution of E_tr and E_s to ET is observed. These years (2008 and 2010) are characterized by the lowest canopy development of our study period (maximum LAI below 4.1 m² m⁻²) and a late growth attributed to the cultivation of a different species.

Similarly, the air temperature that was ~10% higher than that of other years encouraged soil evaporation in 2008 especially in May and June (see Figure 10b). By contrast, 2010 suffered from low temperature values below optimal for quick seed germination after planting. In 2008, 2010, and 2012, soil evaporation was never lower than 10% of ET, even at the vegetation. On the other hand, in September 2014, July 2015, and August 2019, which represent the peak months, E_s consumed less than 3.5% of ET. This behavior is related to the high vegetation fraction, which reduced the available energy for soil evaporation.

At the seasonal scale, 2015 has the highest transpiration of 344 mm and lowest soil evaporation of 136 mm due to dense vegetation favorable temperature and well-distributed rainfall. An examination of the relationship between LAI and E_tr (not shown) revealed a strong positive relationship with correlation coefficient ranging from 0.77 to 0.89 on a daily time scale.

4.3.2. Maize Water Budget

Figure 11 compares the cumulated (P + Irrig), ET, E_tr, and E_S from ISBA-MEB and ΔS.

Table 7. Cumulated water inputs (P + Irrig), measured evapotranspiration (ET_INSITU), simulated transpiration (E_trISBA-MEB), soil evaporation (E_SISBA-MEB), the computed change in the soil water storage (ΔS), and the seepage/extraction (D) beyond 50 cm depth. All measurements are in mm.

Year	2008	2010	2012	2014	2015	2019
P + Irrig	229.1	347.9	343.9	362.7	333.8	368.1
ETINSITU	366.2	439.2	399.7	400.8	477.5	408.5
EtrISBA-MEB	214.1	227.5	288.8	220.4	342.5	302.0
ESISBA-MEB	129.7	180.1	116.8	134.1	60.9	131.9
ΔS	109.2	103.8	68.2	58.3	61.5	136.2
D	−27.9	12.5	11.4	20.2	−82.2	95.

summarizes the water budget components in mm for the study period. The negative and positive sign of D indicates soil water extraction and drainage from depths below 50 cm, respectively.

The striking feature of Figure 11 is that all years behave closely in terms of balance between water inputs and ET. Over all the maize years, (P + Irrig) was lower than ET (Table 7). This behavior has already been observed over several agricultural sites [59,85]. This at first indicates that there is almost no drainage at this site, which is attributed to an impervious soil layer at around 60 cm depth that limits drainage fluxes. Indeed, drainage values are low except in 2019 with a value of 95 mm. This is probably due to the two strong irrigation events in July (events are shown on Figure 5, Figure 6, Figure 7 and Figure 8) scheduled by the farmer to refill the soil reservoir that was severely depleted by the drought period that occurred at this time.

By contrast, the negative drainage value of −82 mm in 2015 could be related to available stagnant water above the impervious layer that contributed to the strong development of the plants during this season (characterized by the highest LAI value (Figure 2) and the highest transpiration amount (Figure 10)). Low drainage values observed in 2010, 2012, and 2014 are characterized by average climatic conditions in terms of evaporative water demand and water inputs. Estimated drainage in 2010, 2012, and 2014 accounted for 12%, 8%, and 12%, while it reached 63% of irrigation amount in 2019. This excess water during the “average” years (2010, 2012, and 2014) is lower than the 22% lost in this region [88].

ΔS values ranging from 58 mm to 136 mm indicate that the plants had extracted a substantial amount of water. This means that a full reservoir is needed at the start of the season to supplement irrigation and rainfall inputs. As no relation is observed between LAI and ΔS in most years, the magnitude of ΔS is controlled by water availability rather than by canopy density. This observation is corroborated by the study of [89], who observed during the growing season of crops, the amount of rainfall influences strongly the amount of soil–water extracted. Another interesting feature is that the cumulative ET values are of comparable magnitude across the years except in 2015 with a value exceeding 470 mm. This is a reconfirmation that the ground water reservoir at the FR-Lam site is key in the estimation of realistic soil water and energy fluxes. It also explains why after each maize season; the site ends up with a very dry soil profile.

5. Conclusions

In the first part, this study compared the performance of a single- and a double-energy budget SVAT model in simulating energy and water fluxes over six seasons of irrigated maize using in situ measurements acquired on a site located in the south west of France. A specific focus was given to the partitioning of the ET into soil evaporation and plant transpiration using sap flow measurements. The ISBA model (single source) and the new multi-energy balance version (ISBA-MEB) were chosen.

• ISBA-MEB provided better estimations of the energy budget components during all the growing seasons, and the main added value of this model was in the prediction of H and to a lesser extent LE during the growth stages when the site is heterogeneous with sparse vegetation. This means that for the future projection of the hydrological functioning and water consumption of maize, a multi-energy balance approach should be preferred to single-source models.
• Concerning the evapotranspiration partition, both ISBA and ISBA-MEB showed good ability to reproduce the observed transpiration derived from sap flow measurements. Nevertheless, the campaign of these measurements took place when the canopy was fully developed, with homogeneous cover.

In the second section, ISBA-MEB was used to analyze the inter-annual variability of the maize water budget as it provided better agreement with the in situ measurements.

• On average, transpiration accounted for a large percentage of ET for most of the years: about 60% for the wet years (2012, 2014, 2015, and 2019) and 46% for the drier years (2008 and 2010). Nevertheless, the partitioning of the ET exhibited a strong year-to-year variability that was closely related to crop development measured by the LAI in this study.
• Another striking feature is that all years behave closely in terms of balance between water inputs and ET; (P + Irrig) was consistently lower than ET except in 2019, a very specific year characterized by a long period of drought in July. This is attributed to very specific conditions on our study site characterized by an existing impervious layer at around 60 cm depth that limits drainage fluxes. This impervious layer, coupled with the good irrigation practice of the farmer, precludes this study from offering recommendations on irrigation scheduling.

As a conclusion, multi-energy approaches should be preferred to single-source models for future projection of water consumption in irrigated maize, even for presumably homogeneous canopies. Nevertheless, further local scale evaluations of the ISBA-MEB model over different crop types and climates are still needed in order to strengthen our conclusions and to feed the global database of land surface parameters such as the ECOCLIMAP II.