Interpreting carbon-water trade-offs in Daisy crop model using Pareto-based calibration

Improving the simulation of carbon and water exchanges is crucial for reliable crop modelling under changing climate conditions. Although model calibration is a key step, optimising multiple outputs can be challenging and often reveals trade-offs between calibration objectives. We applied a Pareto-based multi-objective calibration with the Speed-constrained Multi-objective Particle Swarm Optimisation (SMPSO) algorithm to the Daisy soil–plant–atmosphere model, targeting dry matter (DM), net ecosystem exchange (NEE) and latent heat flux (LE) of winter wheat crops.

The optimal parameter set provided satisfactory results for all objectives (RMSE =0.948 t ha⁻¹ for DM, 1.49 gC m⁻² d⁻¹ for daily NEE and 30.7 W m⁻² for daily LE) but revealed singular trade-offs. The strong compromise between dry matter and NEE likely suggests wrong parameterisation and/or measurement bias, while the trade-off between NEE and LE reflects equifinality issues from evapotranspiration partitioning. Lastly, this analysis also highlighted limitations in simulating stomatal regulation during heatwave conditions, supporting the decoupling between transpiration and carbon assimilation. These findings show that Pareto-based calibration can also serve as a diagnostic tool, identifying structural weaknesses and guiding targeted improvements in process representation for more robust crop model evaluation.

1 Introduction

Crop models are essential tools for yield forecasting and decision-making, but also for assessing the impact of climate change and understanding the complex interactions within the soil–plant–atmosphere continuum (Asseng et al., 2019). Over recent decades, they have been actively developed to account for changing environmental conditions such as rising temperature and elevated atmospheric CO₂ (Timlin et al., 2024). As a result, crop models now integrate more complex approaches and detailed processes, replacing, for instance, the tipping-bucket approach with Richard's equation for soil water dynamics, or using biochemical photosynthesis models (Farquhar et al., 2001) instead of radiation-use efficiency. Biochemical models are often coupled with stomatal conductance models to jointly simulate CO₂ and water fluxes (Keenan et al., 2010). Such couplings are now implemented in advanced crop models such as GECROS (Yin and van Laar, 2005) and Daisy (Plauborg et al., 2010), as well as in most terrestrial biosphere models (e.g., SiB2, ORCHIDEE, IBIS, JULES, etc.). These biosphere models have also been combined with crop models in order to improve crop development and yield predictions while accurately simulating gas exchanges. Some examples are SiBcrop (Lokupitiya et al., 2009), ORCHIDEE–STICS (de Noblet-Ducoudré et al., 2004), Agro–IBIS (Kucharik, 2003) and JULES–SUCROS (Van den Hoof et al., 2011).

In these couplings, carbon and water fluxes between vegetation and atmosphere are assumed to only be regulated by stomatal behaviour. Widely used models of stomatal conductance g_s can be generally expressed as:

$$\begin{array}{}\text{(1)}& g_{\mathrm{s}}=g_{\mathrm{0}}+g_{\mathrm{1}}\cdot f\left(E\right)\cdot A_{\mathrm{net}},\end{array}$$

where A_net is the net leaf CO₂ assimilation rate, f(E) is a function of environmental influences such as Vapour Pressure Deficit (VPD), g₀ and g₁ are parameters. When water is not limiting, plants open their stomata to maximise their carbon uptake. Conversely, they reduce their stomatal aperture under water stressed conditions, compromising between carbon assimilation and water loss (Manzoni et al., 2011). This response is often represented by a decrease of the slope parameter g₁ (Zhou et al., 2013; Liu et al., 2009). However, as stomatal behaviour is influenced by multiple drivers (soil water status, atmospheric dryness, leaf temperature, etc.) which often co-exist during extreme events such as droughts or heatwaves, plant responses may sometimes be contradictory (Grossiord et al., 2020). To date, these complex interactions remain poorly understood and insufficiently represented in models (Liu et al., 2024b; Sabot et al., 2022).

In this context, where stomatal behaviour controls both carbon assimilation and water loss, model calibration requires simultaneous assessment of different outputs. Calibrating against a single output (e.g., evapotranspiration) can result in compensating errors or masking structural deficiencies (Cameron et al., 2022). Multi-objective calibration therefore provides a well-suited framework, although it presents a significant challenge due to the high-dimensional parameter space of such models, the large computational cost of search algorithms and the potentially conflicting nature of targeted objectives. Extensively used in calibration of hydrological models (e.g., Efstratiadis and Koutsoyiannis, 2010; Jahandideh-Tehrani et al., 2020; Moges et al., 2020), Multi-Objective Optimisation (MOO) algorithms aim to optimise all objectives simultaneously and explicitly handle conflicting objectives by exploring trade-offs using a Pareto-based sampling approach (Sharma and Kumar, 2022). This method allows modellers to select any points on the Pareto front, i.e., select an optimal parameter set depending on the modeller's preference (Tang et al., 2018).

Beyond identifying optimal parameter sets, Pareto-based calibration can also provide valuable insights into model behaviour. It can help detect ill-posed models, revealing structural inadequacies, missing processes, or equifinality (i.e., when multiple parameter sets provide similar results; Kollat et al., 2012; Efstratiadis and Koutsoyiannis, 2010). Recently, Harvey et al. (2023) analysed the change of the Pareto front solutions between calibration and validation sets, revealing the influence of data type and model structure. Another study used Pareto-optimal solutions as informative prior distributions to guide posterior uncertainty estimation within a Bayesian framework (Tang et al., 2018).

In this case study, five growing seasons of winter wheat, cultivated on a Belgian site, were simulated using Daisy soil–plant–atmosphere model (Hansen et al., 2012). In addition to the stomatal coupling implemented with Farquhar biochemical model (Plauborg et al., 2010), this open-source model explicitly simulates the surface energy budget. Soil water movement is described by Richards' equation while accounting for drainage pipes, which are widely installed in agricultural fields. Together, these features make Daisy suitable for studying crop yield as well as coupled carbon and water flux dynamics under current and future conditions.

In a previous paper, we identified key parameters controlling crop yield, energy and carbon fluxes and discussed the implications of contrasting water conditions on sensitivity analysis results (Delhez et al., 2025). The objective of this follow-up study was twofold: to improve the simulation of biomass, carbon and water fluxes through calibration, and to explore the trade-offs between these outputs. To this end, we applied a swarm-based MOO algorithm called Speed-constrained Multi-objective Particle Swarm Optimisation (SMPSO) targeting dry matter of vegetation organs (DM), Net CO₂ Ecosystem Exchange (NEE) and latent heat flux (LE). Accurate simulation of these variables is essential for food production assessments and for quantifying crop carbon and water budgets. Furthermore, analysing trade-offs, particularly between carbon (NEE) and water (LE) fluxes, could reveal potential limitations in their coupled representation within the model.

2 Materials and methods

2.1 Case study

The study site is a 12 ha cropland located in Lonzée, Belgium, with a temperate maritime climate ($\stackrel{\mathrm{‾}}{T}=\mathrm{11.1}$ °C and annual precipitation of 716 mm). Cultivated for more than 85 years on a silt loam soil, this site has been equipped with an eddy covariance monitoring system and a micrometeorological station since 2004, as a part of the Fluxnet network. Since 2017, this monitoring site has integrated the ICOS network as a Level 2 station named BE-Lon (https://www.icos-cp.eu/, last access: 5 June 2025).

This study focuses on five growing seasons of winter wheat (Triticum aestivum L.) which are summarised in Table 1. Sown in late October, varieties from S17 and S22 are known to reach heading and maturity stages faster than Sahara and Smart cultivars (Meza et al., 2018, 2023). As the same cultivar was sown during S11 and S15, S11 was selected for the validation. During these five seasons, regular measurements of biomass were performed along with continuous measurements of NEE and LE. These flux data were corrected for frequency losses and filtered for low turbulence and nonstationary conditions, contributing to high-quality coverage detailed in Table 1. Flux data treatment is described in De Ligne et al. (2010) and Pastorello et al. (2020). Data collection (fluxes, biomass and soil) as well as additional information on the site, measurements, or management activities can be found in Delhez et al. (2025).

Table 1 Summary of winter wheat crops in BE-Lon. Coverage of NEE and LE corresponds to the fraction of high-quality hourly data during the calibration period.

2.2 Daisy model

2.2.1 Model description and set-up

Daisy is a crop process-based model simulating carbon, nitrogen and water cycles of the soil–vegetation–atmosphere continuum. Initially published by Hansen et al. (1991), Daisy has since been extensively described (Abrahamsen and Hansen, 2000; Hansen et al., 2012), developed (Plauborg et al., 2010; Holbak et al., 2021), evaluated and compared (Palosuo et al., 2011; Tenreiro et al., 2020; Yin et al., 2020) in the literature. Different modelling options are included in Daisy and must be selected by the modeller depending on the objectives and data availability.

In the present study, the 1D version of Daisy was used with an hourly time resolution. Hourly measurements of air temperature, air vapour pressure, global radiation, wind strength and precipitation were provided as inputs. Water movement in unsaturated soil is described by Richard's equation and, for each soil layer, retention and conductivity curves have been established with the Mualem–van Genuchten model (Mualem, 1976; van Genuchten, 1980). Water drains, located at 70 cm depth in BE-Lon, were explicitly considered (Delhez et al., 2025). Below the soil profile, there is an aquitard layer (i.e., water blocking layer) and an aquifer whose pressure varies according to water table depth measurements also provided as inputs. Potential evapotranspiration is estimated with the FAO56 Penman-Monteith equation, using the reference evapotranspiration (ET0) adjusted to surface cover through bare soil and crop coefficients (Allen et al., 1998; 2006). Water can be evaporated from free water surfaces (canopy interception or ponding), from soil surface or transpired by plants. Actual soil evaporation is governed by the soil hydraulic conductivity, as soil water must be transported to the soil surface, and limited by the potential evaporation. As for the actual transpiration, it is computed together with the entire surface energy balance within the soil–vegetation–atmosphere transfer (SVAT) module (Plauborg et al., 2010). To numerically solve this energy balance, the SVAT module is coupled with the biochemical photosynthesis model of de Pury and Farquhar (1997) and a stomatal conductance empirical model (Eq. 1). In this version, the stomatal conductance g_s is based on Leuning (1995) and adjusted to account for hydraulic signalling (Plauborg et al., 2010):

$$\begin{array}{}\text{(2)}& g_{\mathrm{s}}=g_{\mathrm{0}}+{\displaystyle \frac{g_{\mathrm{1}}\cdot e^{-\mathit{\delta }\left|{\mathit{\psi }}_{\mathrm{c}}\right|}}{\left(\mathrm{1}+\frac{\mathrm{VPD}}{D_{\mathrm{0}}}\right)\cdot \left(C_{\mathrm{s}}-\mathrm{\Gamma }\right)}}{A}_{\mathrm{net}},\end{array}$$

where C_s is the CO₂ concentration at the leaf surface, Γ is the CO₂ compensation point, VPD is the vapour pressure deficit at the leaf surface, ψ_c is the crown water potential (with the related parameter δ) and g₀, g₁ and D₀ are Leuning parameters.

Compared to Daisy v7.0.7, several modifications related to the diffuse radiation, maintenance respiration and carbon reserve remobilisation were made (Appendix A) and later integrated into the latest version (v7.1.0).

2.2.2 Model parameters

In our global sensitivity analysis performed on yield, NEE and LE, we demonstrated the significant influence of about 25 parameters out of 200 (Delhez et al., 2025). However, some adjustments to our selection have been made for this present work. The influential parameters related to the phenology, the allocation and the diffuse fraction of solar radiation have already been estimated separately, as direct measurements are carried out at BE-Lon station. Additionally, the temperature scaling constant c_Vm and the leaf RuBisCO-N fraction were discarded and set to their default value because of their high correlation with other parameters. Conversely, some parameters were added as this study takes new processes into account, such as the photosynthetic contribution of stem and storage organs. Hence, there are 22 parameters included in this study, all described in Table 2. Half of these parameters were considered to be cultivar-specific, while the others are site-specific.

Table 2 Model parameters used in this study and their acceptable range.

^a Parameter estimated separately for each growing season (cultivar-specific). ^b The acceptable ranges were established from measurements, literature and the Daisy reference manual.

2.3 Multi-objective calibration

2.3.1 Key concepts

Calibrating model parameters while considering multiple objectives simultaneously can be formulated as a minimisation problem (Eq. 3):

$$\begin{array}{}\text{(3)}& min\mathit{F}\left(\mathit{x}\right)=\left[\begin{array}{l}{f}_{\mathrm{1}}\left(\mathit{x}\right)\\ f_{\mathrm{2}}\left(\mathit{x}\right)\\ f_{\mathrm{3}}\left(\mathit{x}\right)\end{array}\right],\end{array}$$

where x is a vector of model parameters to be optimised and f_i represents the objective functions. In this study, the objective functions correspond to the relative root mean square error (rRMSE) between modelled and observed values (Eq. 4):

$$\begin{array}{}\text{(4)}& f_i\left(\mathit{x}\right)={\displaystyle \frac{\sqrt{\frac{\mathrm{1}}{n}{\sum }_{k=\mathrm{1}}^n{\left(y_k-d_k\right)}^{\mathrm{2}}}}{\frac{\mathrm{1}}{n}{\sum }_{k=\mathrm{1}}^n{d}_k}},\end{array}$$

Here, y_k and d_k denote the modelled and observed values at time step k for DM (i=1), NEE (i=2) and LE (i=3). Normalising the RMSE by the mean of observed data in Eq. (4) enables fair comparison across variables of different scales.

As these objectives can be conflicting, the minimisation problem defined in Eq. (3) typically results in a set of trade-off solutions known as the Pareto-optimal set, rather than a single solution (Deb, 2001). In this configuration, improving one objective cannot occur without degrading at least one of the others. For this reason, these solutions are also referred to as non-dominated solutions. Visualising these non-dominated solutions in the objective space, referred to as the Pareto front, helps to assess the relationship among objectives and identify model structural errors (Efstratiadis and Koutsoyiannis, 2010).

2.3.2 Search algorithm – SMPSO

Particle Swarm Optimisation (PSO) is a population-based optimisation technique inspired by the social behaviour of birds within a flock (Kennedy and Eberhart, 1995). Each potential solution, called a particle i, is characterised by its position and velocity, and a swarm is defined as the population of particles. Aiming to minimise the objective, the particle moves through the search space, balancing between personal exploration (local optimum) and social learning from the swarm (global optimum; Sharma and Kumar, 2022). The updated position of a particle ${\mathit{x}}_i^{(t+\mathrm{1})}$ is determined by Eq. (5):

$$\begin{array}{}\text{(5)}& {\mathit{x}}_i^{(t+\mathrm{1})}={\mathit{x}}_i^{\left(t\right)}+{\mathit{v}}_i^{(t+\mathrm{1})},\end{array}$$

where ${\mathit{v}}_i^{(t+\mathrm{1})}$ is the updated velocity which explicitly expresses the balance between personal and social learning:

$$\begin{array}{}\text{(6)}& {\mathit{v}}_i^{(t+\mathrm{1})}=w{{\mathit{v}}_i}^{\left(t\right)}+c_{\mathrm{1}}{r}_{\mathrm{1}}\left({\mathit{x}}_{p_i}-{\mathit{x}}_i^{\left(t\right)}\right)+c_{\mathrm{2}}{r}_{\mathrm{2}}\left({\mathit{x}}_{g_i}-{\mathit{x}}_i^{\left(t\right)}\right),\end{array}$$

In Eq. (6), ${\mathit{x}}_{p_i}$ is the best position that particle i has visited, ${\mathit{x}}_{g_i}$ is the best position that the entire swarm has visited (leader), w is the inertia weight controlling the influence of past velocities, r₁ and r₂ are random numbers drawn from the uniform distribution U(0,1), and c₁ and c₂ are learning factors determining the relative influence of personal and global experience.

In this study, we employed the Speed-constrained Multi-objective Particle Swarm Optimisation (SMPSO) algorithm, which is specifically designed for multi-objective problems (Nebro et al., 2009). Based on the OMOPSO algorithm (Sierra and Coello Coello, 2005), SMPSO also uses Pareto dominance (Sect. 2.3.1) and a crowding factor to identify non-dominated solutions. Moreover, it includes a velocity constriction mechanism to prevent the so-called swarm explosion effect, an external archive to store these non-dominated solutions, and the use of polynomial mutation. As suggested by Nebro et al. (2009), the learning factors c₁ and c₂ were drawn from the distribution U(1.5,2.5) and the mutation probability was set as the inverse of the problem's dimensionality.

2.3.3 Application

The Python-based framework for MOO called jMetalPy (Benítez-Hidalgo et al., 2019) includes the SMPSO algorithm and was used in combination with PyDaisy library. For each parameter set (i.e., potential solution) generated by SMPSO, the parameter values were transcribed into Daisy setup files. These include cultivar-specific parameters, with possibly different values for each season, and site-specific parameters that are shared (Table 2). Hence, the model was executed four times, specifically for S15, S17, S19 and S22. The three objective functions (rRMSE) were then computed from these four aggregated runs and passed back to the SMPSO algorithm. Throughout the optimisation process, the hypervolume indicator (referred to as the size of the dominated space) was estimated to evaluate the performance of SMPSO and check for convergence (Guerreiro et al., 2021; Shang et al., 2021). In total, 50 000 model simulations were performed with a swarm size of 200 particles.

2.4 Validation

As mentioned in Sect. 2.1, S11 was not used for this multi-objective calibration but kept for validation. Parameter sets were selected from the Pareto front established with the SMPSO algorithm, and cultivar-specific parameters were taken from S15 as the same cultivar was sown during both seasons. In order to evaluate the performance of Daisy model, four criteria were estimated: (i) the RMSE expressed in the same units as the variable of interest, (ii) the rRMSE for comparison within calibration results, (iii) the Normalised Mean Error (NME) revealing bias and (iv) the model efficiency (EF) which is a normalised metric. This last criterion is also known as Nash–Sutcliffe Efficiency coefficient (NSE) and is computed as follows:

$$\begin{array}{}\text{(7)}& \mathrm{EF}=\mathrm{1}-{\displaystyle \frac{{\sum }_{k=\mathrm{1}}^n{\left(y_k-d_k\right)}^{\mathrm{2}}}{{\sum }_{k=\mathrm{1}}^n{\left(d_k-\stackrel{\mathrm{‾}}{d_k}\right)}^{\mathrm{2}}}},\end{array}$$

In Eq. (7), y_k and d_k denote the modelled and observed values at time step k, and $\stackrel{\mathrm{‾}}{d_k}$ is the mean observed value. The model is perfect when EF = 1, while model predictions are as accurate as the mean of observations when EF = 0, which is easily interpreted. The validation results are considered acceptable if EF is greater than 0.5 (Dumont et al., 2014).

3 Results

3.1 Pareto front

Non-dominated solutions based on the three objective functions are depicted in Fig. 1a. The red point indicates the most balanced compromise, defined as the non-dominated solution nearest to the point representing the lowest values across all objectives. The trade-off between DM and NEE is the most clearly defined, as shown in the 2D projection (Fig. 1b), where the curved Pareto front extends over a large part of the objective space. This 2D front portrays 29 % of the solutions from Fig. 1a, and the other solutions (grey crosses) mainly lie along it. Conversely, the LE–NEE front appears angular and narrower, particularly along the LE axis (Fig. 1c). Combined with the relatively high rRMSE, it highlights the difficulty in accurately simulating observed LE as well as the limited model responsiveness for this objective.

Figure 1 Non-dominated solutions in the objective space: (a) 3D Pareto front according to the rRMSE of DM, LE and NEE, and 2D projection of these solutions considering (b) DM and NEE, (c) LE and NEE. Black dots indicate the Pareto front in 2D space, and the red point highlights the most balanced compromise.

Values of the most balanced parameter set (i.e., red point) can be found in Appendix B (Table B1) along with algorithm convergence indicator (Fig. B1). In Fig. 1a, the specific coordinates of this point are (0.27, 0.40, 0.52) which correspond to RMSE values of 0.948 t ha⁻¹, 3.44 µmol m⁻² s⁻¹ and 62.7 W m⁻² for DM, NEE and LE respectively, computed across all four growing seasons. In the following sections, the simulated outputs using this parameter set are compared to observed values in terms of dry matter as well as carbon and water fluxes.

3.2 Dry matter

In addition to the biomass of individual simulated organs (leaf, stem, storage organs (SOrg) and root), the aboveground vegetative biomass (leaf + stem) was also included in Fig. 2 as distinct measurements of these two organs were not systematic (Table 1). Overall, the simulated dry matter aligns relatively well with observations. Nevertheless, SOrg biomass is consistently underestimated at the end of all growing seasons, resulting in a lower simulated crop yield. Although not shown, predictions of Crop Area Index (CAI) follow a similar pattern, with greater underestimations near the end of the season.

Figure 2 Daily evolution of predicted (solid line) and observed (dot) dry matter for each organ (t ha⁻¹). Error bars depicted for the observed dry matter represent the standard deviations and the dashed vertical line indicates the anthesis.

For S15, this underestimation may be attributed to insufficient remobilisation of stem reserves, as indicated by the high simulated stem biomass in July, even greater than the combined leaf and stem observations. However, this explanation does not hold for the other seasons, where stem biomass appears realistic. Beyond the capacity of remobilised reserves, parameters such as the efficiency of this remobilisation (considered fixed in Daisy) as well as contributions from stem and ear to photosynthesis (SOrgPhotEff and StemPhotEff) may also affect SOrg biomass during the reproductive stage (after anthesis).

3.3 Fluxes

As depicted in Fig. 3, simulated temporal variations follow observed trends, although predictions tend to be more negative (i.e., more CO₂ uptake) during the vegetative stage. This underestimation is particularly evident for S17, also highlighted by the scatter plot, where the regression slope is significantly higher than 1 (β=1.14). This bias suggests an overestimation of photosynthesis and/or an underestimation of respiration. Despite this discrepancy, the model accurately simulates the biomass of vegetative organs during this stage. On the other hand, during the reproductive stage (when SOrg is underestimated; Fig. 2), NEE predictions no longer show a clear bias. This partly explains the observed trade-off between DM and NEE (Fig. 1b). For instance, increasing the remobilisation of stem reserves during S15 would likely result in a more accurate crop yield, but also in an undesirable increase in NEE as CO₂ is released through this process.

Figure 3 Daily NEE (gC m⁻² d⁻¹) for each growing season (left) with the red shaded area representing standard deviations of observations. Hourly predictions against observed values of NEE (µmol m⁻² s⁻¹; right) where observed values are displayed on the x-axis. The p-value indicates if the regression slope β (red) is significantly different from the 1:1 line (black) and ρ is the Pearson correlation coefficient.

It is also worth noting that Daisy simulates a sudden decrease in NEE at the end of each season, which is not observed. When reaching maturity (DS = 2), plant processes are completely stopped, so is the autotrophic respiration. Actually, matured plants continue to respire and contribute to carbon fluxes until harvest.

As expected from the rRMSE values observed in Fig. 1, LE fluxes are not well captured by the model. Observed daily fluxes exhibit greater variations, frequently reaching high values or dropping below 50 W m⁻² while simulated values remain relatively stable (Fig. 4). This is also illustrated in the scatter plots, where hourly observed data easily reach 600 W m⁻², but simulated fluxes never exceed 500 W m⁻², resulting in a poor linear agreement between observed and simulated fluxes as indicated by the low value of Pearson correlation coefficients ρ (Fig. 4). This disagreement is particularly visible in mid-April 2019, with an observed daily average around 150 W m⁻² in comparison to a prediction of only ∼100 W m⁻². During this period, simulated evapotranspiration equals potential evapotranspiration, suggesting that the plant was not water-limited and that the underestimation stems from deficiencies in energy partitioning rather than water availability.

Figure 4 Daily average of LE (W m⁻²) for each growing season (left) with the red shaded area representing standard deviations of observations. Hourly predictions against observed values of LE (right) where observed values are displayed on the x-axis. The p-value indicates if the regression slope β (red) is significantly different from the 1:1 line (black) and ρ is the Pearson correlation coefficient.

Even though there was no apparent soil water deficit, this period (18–22 April) was still characterised by elevated temperatures and a higher vapour pressure deficit (VPD) compared to preceding days. Weather data during this period are depicted in the Appendix C (Fig. C1). It does not meet the formal Belgian definition of a heatwave, but the daily mean temperatures were more than 5° C above the 1981–2010 averages. During this 5 d period, the observed stomatal conductance g_s reached extremely high values (Fig. 5a), associated with an important increase in LE fluxes (Fig. 5c), while the magnitude of the carbon uptake (GPP) was not affected (Fig. 5b). These observations reveal a decrease in water use efficiency during this period, indicating that transpiration was intensified with no corresponding increase in carbon assimilation. This trend was not captured by the model, as simulated g_s remained lower than 0.03 m s⁻¹ (Fig. 5a). The slight increase in predicted g_s during these five days can be explained by the minor increase in carbon uptake (Fig. 5b).

Figure 5 Evolution of predicted and observed (a) stomatal conductance g_s, (b) Gross Primary Productivity (GPP), (c) LE fluxes and (d) Soil Water Content (SWC) at three different depths. Solid lines are predictions. Observed g_s was deduced from Penman–Monteith equation, measured fluxes and weather data (Appendix C). The axis break indicates different scales but does not imply omitted observations.

In addition to the insufficient simulated transpiration, soil evaporation may also be underestimated. During these five days, soil water content at the surface (5 cm depth) decreased by 6.2 % whereas the model predicted a decrease of 4.2 % (Fig. 5d). We assume that this difference (2 %) is mainly due to an underestimation of soil evaporation rather than drainage as it is limited to the surface layer and because of the particular meteorological conditions. Considering a surface layer from 0 to 10 cm, this is equivalent to a daily contribution of 10.4 W m⁻². Hence, the underestimation of soil evaporation can partly explain the difference in LE fluxes observed in Fig. 4.

3.4 Validation

When using the most-balanced parameter set for S11, NEE is well reproduced by the model with EF greater than 0.8, rRMSE around 40 % and NME close to 0 for both time scales (Table 3a). Although rRMSE and EF also present acceptable results for DM, it is worth noting that there is an overall underestimation of this variable (NME $=-\mathrm{0.222}$).

Table 3 Evaluation criteria for S11 using (a) the most-balanced parameter set (red point) and (b) the parameter set minimising F_LE. RMSE are expressed in W m⁻² for LE, in µmol m⁻² s⁻¹ for hourly NEE, in gC m⁻² d⁻¹ for daily NEE and in t ha⁻¹ for dry matter (DM). Other criteria have no units.

On the other hand, there is no good agreement between observed and predicted LE, especially when considering a daily time scale, as shown by the four indicators. Even when using the parameter set that minimises F_LE (Fig. 1a), LE remains poorly captured by the model (Table 3b). As expected from the Pareto front, this slight improvement in LE predictions is accompanied by a degradation of NEE and DM predictions.

4 Discussion

4.1 Pareto-based calibration

Considering all objectives, the chosen parameter set (red point in Fig. 1) shows satisfactory calibration results, with a total RMSE of 0.948 t ha⁻¹, 3.44 µmol m⁻² s⁻¹ and 62.7 W m⁻² for DM, NEE and LE, respectively. Compared with other models, calibrated RMSE values for total biomass are typically between 0.5 and 1.5 t ha⁻¹ (Zhang et al., 2020; Dumont et al., 2014; Casanova and Judge, 2008). Regarding NEE, reported RMSE values for daily fluxes range from 1.3 to 1.9 gC m⁻² d⁻¹ (Revill et al., 2019; Senapati et al., 2018; Lehuger et al., 2010). In the present study, the RMSE computed from hourly NEE (3.44 µmol m⁻² s⁻¹) corresponds to a RMSE of 1.49 gC m⁻² d⁻¹ when using daily fluxes, which falls in the mentioned range. Casanova and Judge (2008) reported a RMSE of 62.45 W m⁻² for hourly LE fluxes, and other studies presented values ranging from 24.65 to 48.5 W m⁻² for daily fluxes (Wang et al., 2024; Liu et al., 2020a; Dutta et al., 2016), while we found 30.74 W m⁻² using daily averages.

The calibrated parameters are all within the explored range, without relying on extreme values (Table B1), indicating a physically plausible solution. It is, however, worth noting that some cultivar-specific parameters have very different values between seasons, such as δ describing the hydraulic signalling for stomatal conductance or the photosynthetic efficiency of non-leaf organs (SOrgPhotEff and StemPhotEff). Such differences may reflect genotypic variability, but they may also be due to environmental factors not represented in the model. This is, for example, the case for the remobilisation reserve ShldResC which can be impacted by genotype and water stress (Ehdaie et al., 2006). The absence or misrepresentation of environmental factors in the model can be compensated by wider variations between cultivars. Depending on the modeller's objective, these parameters might be kept as cultivar-specific while investigating the environmental influences, or constrained across cultivars to improve model parsimony and robustness.

In addition to providing reliable solutions, Pareto-based multi-objective calibration also offers insights into internal model limitations (Kollat et al., 2012). As related by Schoups et al. (2005), the shape of the Pareto front is an indicator of model structural error. Angular fronts (e.g., NEE–LE front) suggest that different objectives can be simultaneously optimised whereas significant trade-offs (curved or linear fronts) may indicate a wrong parameterisation of the model (Wöhling et al., 2013; Efstratiadis and Koutsoyiannis, 2010). The following sections focus on the NEE–DM and NEE–LE trade-offs and discuss whether the disagreements between observations and model predictions come from uncertainty due to measurement errors or deeper structural issues.

4.2 Where does the carbon go?

As shown by the major trade-off in Fig. 1b, adjusting model parameters to fit NEE fluxes inevitably compromises the simulation of organ DM. The analysis in Sect. 3.2 and 3.3 reveals that this would specifically result in an underestimation of biomass. DM during the vegetative stage is well reproduced when simulated NEE is overly negative, and SOrg DM is underestimated when NEE predictions show less bias (Figs. 2 and 3). This was also observed in the validation results, with a consistent underestimation of DM while NEE showed neither over- nor underestimation. Several factors may contribute to this mismatch:

• Underestimation of heterotrophic respiration. If so, as respiration is regarded as positive and photosynthesis as negative with the micro-meteorological sign convention, NEE predictions would be too negative (as observed). Improving heterotrophic respiration would affect NEE predictions but its influence on DM would be limited, which could weaken the trade-off between these variables. At this study site, Suleau et al. (2011) partitioned soil respiration into its heterotrophic and autotrophic components using soil chambers, from which they derived temperature functions. When comparing their fitted regression to simulated heterotrophic respiration, Daisy appears to underestimate both temperature sensitivity and baseline respiration rate for three growing seasons (Fig. D1), supporting this first assumption. It should be noted, however, that this regression depends solely on temperature, ignoring the seasonal variation in residue and fertilisation inputs that also affect soil respiration dynamics. Furthermore, heterotrophic respiration is most likely to be underestimated as, during the bare soil period before sowing, NEE (which only consists of soil respiration at that time) is mostly underestimated (Table D2). Adjusting soil parameters might reduce this trade-off, but it would require long-term simulations and respiration data from this study site.

• Wrong modelling assumptions. Modelling real-world systems comes with simplifications. For instance, in this study, the partitioning schemes of assimilates were derived from S22, since root measurements were only made during that season, as well as early distinction between leaf and stem. These partitioning coefficients were also applied to the other cultivars, but they might be inadequate and partly contribute to this trade-off. If too much assimilated carbon is allocated to the roots at the expense of the stem, the model would underestimate stem DM, which is crucial for grain filling during the reproductive stage. Conversely, as leaves are not directly affected and are the major contributor to photosynthesis, the influence on NEE would be attenuated. Besides this example, many other assumptions could play a role in these results, such as the organ death rates depending only on development stages or the non-dynamic leaf abscission rate. However, accounting for these within the model would be, in most cases, ineffective as model users usually do not have enough data to support this level of complexity.

• Uncertainties and errors in NEE observations. Eddy covariance (EC) data are affected by multiple factors, including instrumental noise, frequency losses and turbulence conditions (Aubinet et al., 2012). To select only high-quality data, they are thoroughly filtered, resulting in low data coverage, as in Table 1. Moreover, turbulence conditions are accounted for in the uncertainty estimation rather than affecting the data itself, with an uncertainty related to the friction velocity, an indicator of turbulence strength (Pastorello et al., 2020). However, an additional source of errors has been discussed in relation to the energy balance closure problem (Gao et al., 2019). EC systems show a systematic underestimation of energy fluxes (sensible and latent heat) ranging from 15 % to 20 %, which is generally attributed to sub-mesoscale organised structures, generated in the daytime convective layer and enhanced by underlying surface heterogeneities (Mauder et al., 2020; Aubinet et al., 2012). While a standard correction method exists for energy fluxes (Pastorello et al., 2020), its application to CO₂ fluxes remains under active research (Mauder et al., 2024). Some studies suggest that CO₂ fluxes may also be underestimated (Gao et al., 2019; Mauder et al., 2010), though Liu et al. (2024a) found that the bias rather depends on soil water conditions. At BE-Lon, this potential underestimation of NEE is supported by the crop carbon budget previously estimated (Buysse et al., 2017). When combining NEE and field measurements (harvest, residues and inputs), they found that the study site continuously loses an unexpectedly large amount of carbon, depleting the soil much faster than recently observed (Dumont et al., 2025). This highlights that, beyond random error estimates, systematic biases in NEE measurements may propagate into model calibration and interpretation. As Beven (2019) emphasises, observational datasets should not be treated as absolute truth, and acknowledging these uncertainties is important when designing calibration frameworks or interpreting model evaluations.

4.3 Structural limitations in simulating LE

When there is no evaporation from free water surfaces, the surface water balance in Daisy can be summarised as:

1. computation of potential evapotranspiration (PET) using the FAO56 equation

2. partitioning of PET between potential soil evaporation and plant transpiration according to canopy cover

3. computation of actual soil evaporation, depending on potential soil evaporation and soil water availability

4. numerical resolution of Farquhar–Leuning models simultaneously with the SVAT module (surface energy balance) to obtain actual transpiration, surface temperatures and carbon assimilation.

Model limitations may arise from this workflow in two ways. First, the evapotranspiration (ET) is constrained by FAO56 PET which assumes a fixed stomatal conductance (Allen et al., 2006). Although it is widely used in crop modelling, this simplification can reduce the accuracy of ET estimations (Liu et al., 2023; Ghiat et al., 2021). By considering a constant g_s for PET, simulated ET can be severely limited when stomata are actually wide open. Conversely, when stomata are closed, PET is overestimated and can therefore lead to higher soil evaporation than observed. This may partly explain the lack of temporal variation observed in Fig. 4 and the high values of rRMSE (Fig. 1c).

Second, soil evaporation (E_soil) is computed outside the iteration loop that couples the SVAT module with the Farquhar–Leuning models. It is constrained either by soil water availability or by the potential evaporative demand, depending on which is limiting. By contrast, fully coupled models such as the CLM model calculate resistance-based E_soil simultaneously with transpiration and sensible heat fluxes (Swenson and Lawrence, 2014). This ensures that E_soil is driven by microclimate conditions (e.g., canopy air vapour pressure) and reinforces the mechanistic coupling between the surface energy and water balances. However, the appropriate formulation of resistance-based E_soil remains an active research area, as recent studies continue to refine its parameterisation and physical basis (Schulz and Vogel, 2020; Lehmann et al., 2018).

Beside these workflow aspects, the relatively narrow Pareto shape (Fig. 1c) suggests the presence of equifinality, i.e., multiple parameter sets produce equally good results (Her and Seong, 2018). This likely reflects compensatory interactions between soil evaporation and plant transpiration. Since PET is partitioned by canopy cover, parameter changes that alter canopy development can redistribute PET between the two components while leaving total PET unaffected. Under non-stressed water conditions, the ecosystem would meet this evaporative demand, and total ET (i.e., LE outputs) would remain unchanged as well. This masking effect underscores the well-known challenge of ET partitioning in ecosystem modelling, further complicated as eddy covariance systems do not distinguish between soil evaporation and transpiration (Berg and Sheffield, 2019; Stoy et al., 2019; Scott and Biederman, 2017). Recent methods, such as flux mapping, have been developed to better understand model internal fluxes leading to equifinality (Khatami et al., 2019).

As depicted by Fig. 5, discrepancies in LE predictions were especially evident under elevated temperatures and high VPD. While these may partly stem from limitations mentioned above, they also point to deeper issues in the coupling between Farquhar and Leuning models and particularly in how it responds to atmospheric drivers.

4.4 Coupling of carbon and water exchanges

Stomatal behaviour plays a central role in plant regulation of water loss and carbon uptake, and has therefore been extensively studied (Damour et al., 2010). While it is well established that mild edaphic drought induces a decrease in carbon assimilation by reducing the stomatal aperture (Beauclaire et al., 2024), plant responses to heatwaves and high atmospheric demand are less evident. Most studies report reduced g_s under high VPD (Grossiord et al., 2020; Bourbia and Brodribb, 2024), aligning with the theoretical predictions of stomatal models (Sabot et al., 2022). On the other hand, high temperatures seem to trigger stomatal opening to promote transpiration cooling (Urban et al., 2017), but reported results vary greatly among species (Wang et al., 2026; Moore et al., 2021). Under the combined effect of high temperature and VPD, Marchin et al. (2022) observed a significant increase in g_s for two well-watered species, interpreting this as a strategy to prevent leaf overheating via enhanced LE loss without any increase in carbon assimilation. They later confirmed that this behaviour, referred to as stomatal decoupling, is not restricted to well-watered species (Marchin et al., 2023).

This stomatal decoupling aligns with ecophysiological studies, highlighting the role of non-stomatal factors during photosynthesis. Inside the leaf, CO₂ diffuses through air spaces and membranes to reach the sites of carboxylation inside chloroplasts, whereas H₂O moves from the xylem network to the stomata (Flexas et al., 2012; Sack and Holbrook, 2006). These different pathways can be affected independently, where carbon assimilation can be restricted by non-stomatal factors, often dominated by mesophyll conductance (Flexas et al., 2012). Over the last decade, these non-stomatal limitations have received increased attention, especially under water-stressed conditions (Nadal and Flexas, 2018; Gago et al., 2020). Based on leaf gas exchange measurements, experimental studies proved that these limitations should be considered in leaf models linking photosynthesis and transpiration (Yang et al., 2019; Beauclaire et al., 2024; Chen et al., 2025). Despite that, these complex processes are still not fully understood, and hence barely integrated into crop or terrestrial models (Vidale et al., 2021). As an example, an empirical equation for the mesophyll conductance was implemented in the land surface model JSBACH (Knauer et al., 2019). The authors extensively discussed the environmental drivers, the modelling implications, as well as the effects on the predicted fluxes in a subsequent paper (Knauer et al., 2020).

These non-stomatal limitations, or stomatal decoupling, are not modelled in Daisy. Based on stomatal coupling (Eq. 2), the increase in VPD reduced the predicted slope between carbon assimilation and g_s. Consequently, the model tended to underestimate g_s and LE fluxes and overestimated sensible heat, leading to a misrepresentation of energy partitioning under such conditions. On the other hand, observations show that plants likely opened their stomata, promoting water loss, while carbon assimilation did not increase accordingly. This stomatal decoupling was particularly visible in April 2019, as the heatwave lasted several days, but also appeared during shorter warm periods, suggesting a recurrent plant strategy. Future work should (i) continue the effort in better understanding the mechanisms of non-stomatal limitations affecting carbon assimilation and (ii) incorporate these mechanisms in process-based models, particularly under such conditions.

5 Conclusion

In this study, we employed a Pareto-based multi-objective calibration approach to optimise the Daisy soil–plant–atmosphere model, focusing on three key objectives (DM, NEE and LE). The SMPSO algorithm successfully identified parameter sets that balanced these objectives within realistic bounds. Furthermore, this methodology also offered valuable insights into model behaviour and inherent limitations, as particular trade-offs between the objectives were revealed. Notably, the strong compromise between DM and NEE suggests possible wrong parameterisations of the model and/or biases in observational EC data. Moreover, the LE–NEE trade-off indicates issues of equifinality, which underlines the challenging partitioning of evapotranspiration. The inaccuracy of LE predictions was also discussed, underscoring the oversimplification of the FAO56 Penman-Monteith method.

Additionally, the model struggled to reproduce observed increases in transpiration during heatwave conditions, highlighting deficiencies in simulating stomatal responses and in the coupling between carbon and water fluxes. These limitations are particularly relevant given the increasing frequency and intensity of heatwaves under climate change, which can significantly impact ecosystem dynamics. The influence of heatwaves and underlying mechanisms of stomatal decoupling need to be further investigated to improve the simulation of carbon and energy fluxes.

Overall, this work illustrates how multi-objective calibration can go beyond parameter estimation to critically assess model assumptions and guide future improvements. Addressing the structural limitations identified here – particularly in ET partitioning and (non-)stomatal regulation – could strengthen the ability of Daisy and similar crop models to capture ecosystem responses under changing climate conditions.

Appendix A: Daisy modifications

• Diffuse fraction. The diffuse fraction of global radiation f_d was estimated using the optical air mass and assuming a cloudless sky (Appendix 1 of de Pury and Farquhar, 1997). A second modelling option for estimating f_d, adapted for cloudy skies, was implemented. Based on Ridley et al. (2010), f_d depends on the hourly and daily clearness index (resp. k_t and K_t), the apparent solar time (AST) and the solar angle α:

  $$\begin{array}{}\text{(A1)}& f_{\mathrm{d}}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+e^{{\mathit{\beta }}_{\mathrm{0}}+{\mathit{\beta }}_{\mathrm{1}}{k}_t+{\mathit{\beta }}_{\mathrm{2}}\mathrm{AST}+{\mathit{\beta }}_{\mathrm{3}}\mathit{\alpha }+{\mathit{\beta }}_{\mathrm{4}}{K}_t}}},\end{array}$$

  As both diffuse and global radiation are measured at BE-Lon, the five β parameters have been estimated with the Nelder–Mead method, considering the RMSE metric as an objective function. This method has been applied to 10 years of measurements (2014–2023), covering the entire duration of the case study. The resulting parameters (β₀, β₁, β₂, β₃, β₄) are respectively equal to −4.8094, 5.4735, 0.020702, 0.099970 and 1.3944 (RMSE for f_d= 0.164).

• Maintenance respiration. Respiration required to maintain existing biomass is defined as maintenance respiration. In Daisy, the maintenance respiration of an organ $R_{\mathrm{m}}^{\mathrm{org}}$ is proportional to its dry weight W^org:

  $$\begin{array}{}\text{(A2)}& R_{\mathrm{m}}^{\mathrm{org}}=r_{\mathrm{m}}^{\mathrm{org}}\cdot f\left(T\right)\cdot W^{\mathrm{org}},\end{array}$$

  where $r_{\mathrm{m}}^{\mathrm{org}}$ is the maintenance respiration coefficient of the organ and f(T) is a temperature function. In addition to temperature dependence, it was demonstrated that mature or senescing tissues required lower maintenance respiration than younger tissues (Boote et al., 2013). Maintenance respiration declines near the end of the growing season due to a significant fraction of dead cells and metabolically inactive long-term storage (Amthor, 2025). To account for this, a second dependency f(DS) was added into Eq. (A2) and is a user-defined function of development stage. This more realistic dependence was added to verify whether its consideration impacted the trade-off between DM and NEE, which was not the case.

• Remobilisation. After anthesis, carbon reserves stored in the stem can be mobilised to support grain filling. The remobilisation equations in Daisy were modified to account for both remobilisation efficiency and the temporal dynamics of the remobilisation rate. This rate peaks around 20 d after anthesis before declining, resulting in a sigmoid-shaped depletion of stem reserves (Liu et al., 2020b; Ehdaie et al., 2008). Similarly, this process was also modified to better match observations and verify whether it impacted the trade-off between DM and NEE, which was not the case.

Appendix B: SMPSO results

Table B1 Most balanced parameter set among non-dominated solutions (red point in Fig. 1).

Figure B1 Hypervolume indicator with increasing iterations. The red line corresponds to the solutions presented in this study and the grey line is a previous attempt with less iterations (200 iterations, hence 40 000 simulations).

Appendix C: Determination of stomatal conductance

According to the Penman–Monteith equation, the latent heat flux (LE) can be expressed as follows:

$$\begin{array}{}\text{(C1)}& LE={\displaystyle \frac{\mathrm{\Delta }\left(R_{\mathrm{n}}-G\right)+\frac{{\mathit{\rho }}_{\mathrm{a}}{c}_{\mathrm{p}}\mathrm{VPD}}{r_{\mathrm{a}}}}{\mathrm{\Delta }+\mathit{\gamma }\left(\mathrm{1}+\frac{r_{\mathrm{s}}}{r_{\mathrm{a}}}\right)}},\end{array}$$

where R_n is the net irradiance (W m⁻²), G is the ground heat flux (W m⁻²), VPD is the Vapour Pressure Deficit (Pa), ρ_a is the dry air density (equal to 1.22 kg m⁻³), c_p is the specific heat capacity of air (equal to 1013 J kg⁻¹ °C⁻¹) and r_s is the stomatal resistance ($r_s=g_{\mathrm{s}}^{-\mathrm{1}})$. The slope of the saturation vapour pressure with temperature Δ can be estimated based on vapour pressure at saturation e^∗(T_air):

$$\begin{array}{}\text{(C2)}& \mathrm{\Delta }={\displaystyle \frac{\mathrm{4098}{e}^{\ast }\left(T_{\mathrm{air}}\right)}{{\left(T_{\mathrm{air}}+\mathrm{237.3}\right)}^{\mathrm{2}}}},\end{array}$$

Figure C1 Hourly evolution of latent heat flux predictions (grey) and observations (red) along with meteorological data (incoming shortwave radiation, air temperature and VPD).

The psychometric constant γ is a function of atmospheric pressure P:

$$\begin{array}{}\text{(C3)}& \mathit{\gamma }={\displaystyle \frac{c_{\mathrm{p}}P}{\mathrm{0.622}\cdot \mathrm{2.45}\cdot {\mathrm{10}}^{\mathrm{6}}}},\end{array}$$

And finally, the atmospheric resistance r_a is composed of the aerodynamic and boundary layer resistances:

$$\begin{array}{}\text{(C4)}& r_{\mathrm{a}}={\left({\displaystyle \frac{u}{u_{\ast }}}\right)}^{\mathrm{2}}+\mathrm{6.2}{u}_{\ast }^{-\mathrm{0.67}},\end{array}$$

As LE, G, R_n, VPD, the air temperature T_air, P, the wind velocity u, the friction velocity u_∗ are directly measured or at least computed within the ONEFlux pipeline (Pastorello et al., 2020), the stomatal resistance can be determined by combining Eqs. (C1)–(C4).

Appendix D: Heterotrophic respiration

Figure D1 Daily evolution of heterotrophic respiration (gC m⁻² d⁻¹) – model predictions (grey) and estimations based on the temperature function fitted by Suleau et al. (2011; red).

Table D1 Evaluation criteria regarding daily NEE during the bare soil period before sowing. RMSE and ME are expressed in gC m⁻² d⁻¹ and other criteria have no units. S11 was simulated with the most-balanced parameter set, and cultivar-specific parameters were taken from S15.