Development of the DO₃SE-Crop model to assess ozone effects on crop phenology, biomass, and yield

A substantial body of empirical evidence exists to suggest that elevated O₃ levels are causing significant impacts on wheat yields at sites representative of highly productive arable regions around the world. Here we extend the DO₃SE model (designed to estimate total and stomatal O₃ deposition for risk assessment) to incorporate a coupled A_net–g_sto model to estimate O₃ uptake; an O₃ damage module (that impacts instantaneous A_net and the timing and rate of senescence); and a crop phenology, carbon allocation, and growth model based on the JULES-crop model. The model structure allows scaling from the leaf to the canopy to allow for multiple leaf populations and canopy layers. The DO₃SE-Crop model is calibrated and parameterised using O₃ fumigation data from Xiaoji, China, for the year 2008 and for an O₃-tolerant and sensitive cultivar. The calibrated model was tested on data for different years (2007 and 2009) and for two additional cultivars and was found to simulate key physiological variables, crop development, and yield with a good level of accuracy. The DO₃SE-Crop model simulated the phenological stages of crop development under ambient and elevated O₃ treatments for the test datasets with an R² of 0.95 and an RMSE of 2.5 d. The DO₃SE-Crop model was also able to simulate O₃-induced yield losses of ∼11 %–19 % compared to observed yield losses of 12 %–34 %, with an R² of 0.68 (n=20) and an RMSE of 76 g m⁻². Additionally, our results indicate that the variance in yield reduction is primarily attributed to the premature decrease in carbon assimilation to the grains caused by accelerated leaf senescence, which is brought forward by 3–5 d under elevated O₃ treatments.

1 Introduction

Ground-level ozone (O₃) is considered the most critical air pollutant causing global damage to agricultural crops. Elevated O₃ concentrations are particularly problematic in Asia, where decades of rapid economic growth, industrialisation, and urbanisation have seen sharp rises in pollutant emissions associated with burning fossil fuels (Lin et al., 2017), causing substantial O₃-induced crop yield losses across the region (Z. Feng et al., 2022). At the same time, climate change is considered a substantial threat to arable productivity through changes in average and extreme temperature and precipitation profiles across the region (IPCC, 2021). Reductions in precipitation are considered responsible for poor harvests in recent years (Liu et al., 2010), and rising temperatures that reduce the length of the crop growing season are thought to have caused losses in crop yield (Malhi et al., 2021). There is now substantial evidence showing that stresses from O₃ pollution and climate variability interact, causing either additive, synergistic, or antagonistic responses in crop development, growth, and yield (Sillmann et al., 2021). The threat posed by these stresses is a particular cause for concern in Asia since the continent contributes approximately 43 % of the global wheat production, with China contributing the highest production levels at 17 % of the global wheat supply (Feng et al., 2021). O₃ levels are rising substantially in important wheat-growing areas in China such as the North China Plain and the Yangtze River Delta (Li et al., 2020; Zhang et al., 2023). Concern over O₃ impacts led to the implementation in 2013 of a range of policies to try to reduce O₃ precursor emissions across China. These included a comprehensive management plan to control volatile organic compounds (VOCs) from key industries, an atmospheric pollution prevention and control law of the People's Republic of China, and a 2020 VOC management plan (Li et al., 2021). As a result, nitrogen oxide (NO_x) emissions, an important O₃ precursor, have decreased by 21 % from 2013 to 2017 (Li et al., 2021). By contrast, VOCs have only slightly decreased by 2 % over the same period. Since China has a VOC-limited O₃ regime, the reductions in NO_x lead to rather insignificant changes in O₃ concentration (Li et al., 2021), though evidence suggests that reductions in O₃ may be higher in rural areas than in urban areas (Lee et al., 2020). This implies future policies to tackle ground-level O₃ pollution in China need to increase their focus on reducing VOCs along with NO_x (Lee et al., 2020) and also emphasise the importance of being able to make assessments of O₃ damage to key receptors such as staple crops.

At present, methods to assess the risk to crop productivity from changes in O₃ and climate variables use a variety of different O₃ risk assessment methods (Ronan et al., 2020) and crop models as discussed in depth in Emberson et al. (2018). In the past, O₃ risk assessment methods relied heavily on dose–response relationships, empirically derived relationships that assess changes in a response variable (most commonly yield) against an O₃ exposure metric (concentration or, more recently, flux-based indices) (Pleijel et al., 2022). By contrast, methods to assess the impact of climate variables (most commonly changes in temperature, precipitation, and CO₂ concentration) tend to use crop models since these allow the integration of the combined effect of a number of different variables acting simultaneously to affect crop development, growth, and yield (Schauberger et al., 2019). A new generation of crop models that include O₃ damage are now being developed and applied and have the potential to estimate the combined effect of O₃ and climate variables on crop development, biomass, and yield. Such models can arguably be classified into two types of crop model: firstly, those that rely on O₃ metrics (e.g. AOT40 or M7) to modify crop growth determined by radiation use efficiency (Guarin et al., 2019, 2024) or evapotranspiration (Droutsas et al., 2020); secondly, those that estimate stomatal O₃ uptake to modify crop growth determined by photosynthesis and subsequent carbon assimilation (Tao et al., 2017; Schauberger et al., 2019; Nguyen et al., 2024). The DO₃SE-Crop model falls into the latter category of photosynthesis-based crop models and was developed to bridge the gap between O₃ risk assessment modelling methods and crop models.

The DO₃SE model is an O₃ deposition model that can be embedded within atmospheric chemistry transport models (e.g Simpson et al., 2012) and uses either a multiplicative or coupled A_net–g_sto model to estimate stomatal O₃ flux (Pande et al., 2024). The accumulated stomatal O₃ flux has been successfully used as a damage metric (POD_y – phytotoxic ozone dose over a threshold y; LRTAP, 2017) to predict O₃-induced yield loss (Pande et al., 2024). The ability of the DO₃SE model to simulate A_net, as well as the inclusion of a process-based O₃ damage module for both instantaneous A_net and early and enhanced senescence (after Ewert and Porter, 2000), lends itself to the development of the DO₃SE model as a process-based crop model. The inclusion of resistance algorithms that can assess the transport of O₃ concentrations from a reference height above a canopy down to the canopy top means the model can be embedded within existing atmospheric chemistry transport schemes and hence applied to regional- or global-scale O₃ risk assessment whilst also modelling O₃ deposition. A comparison of the coupled stomatal A_net–g_sto model with the multiplicative g_sto model within the DO₃SE framework was made in Pande et al. (2024), and it showed that the A_net–g_sto model performed equally well, if not better, when used to develop O₃ dose–response relationships for European wheat. This provides evidence of the suitability of the new photosynthesis-based g_sto model in DO₃SE.

In this study, we describe the development of a new DO₃SE-Crop model which builds on the modified stomatal deposition component of the DO₃SE model (Pande et al., 2024) so that both CO₂ uptake for carbon assimilation and O₃ uptake via the stomata can be modelled consistently. Further, we have incorporated the UK JULES-crop model (Osborne et al., 2015) to allocate assimilated carbon to plant components (roots, leaves, stems, and harvest organs) according to crop development stage. We also take account of the modifying effect of O₃ on instantaneous A_net as well as accumulated A_net via O₃ effects on the onset and rate of leaf senescence and timing of crop maturity through incorporation of algorithms developed by Ewert and Porter (2000). The UK JULES-crop model is used since this is the UK land surface exchange scheme in the UK Earth System Model (UKESM) (Osborne et al., 2015) which has recently been developed to include exchange and impact of trace gases (including O₃) along with other biogeochemical cycling between the atmosphere and the land surface (Leung et al., 2020). This would in the future allow comparison of the UK JULES-crop model, which uses O₃ mechanisms that modify instantaneous A_net to mimic changes in yield consistent with flux–response relationships (Sitch et al., 2007), with the alternative O₃ damage mechanisms used within DO₃SE-Crop.

Figure 1 Schematic of the DO₃SE-Crop model.

Here, we calibrate and evaluate the DO₃SE-Crop model using an experimental FACE dataset collected in Xiaoji, China. This allows us to investigate the ability of the model to simulate O₃ damage for a comparable agro-ecological region where crop productivity is severely threatened by both O₃ pollution and climate change. The key objectives of the paper are to assess the ability of DO₃SE-Crop to simulate (i) key phenological stages, (ii) the relationship between leaf-level physiological variables and within-canopy O₃ concentrations, (iii) C allocation to different parts of the crop, and (iv) O₃-induced yield losses for tolerant and sensitive cultivars.

2 Methods

2.1 DO₃SE-Crop model

Here we describe the development of the DO₃SE-Crop model. In this study, version 4.39.16 of the DO₃SE-Crop model was used (Bland, 2024) for wheat (Triticum aestivum), which is widely considered to be one of the most sensitive staple crops to O₃ (Feng et al., 2018). The key components of DO₃SE-Crop are illustrated in Fig. 1. The model integrates meteorological data, crop parameters, and site characteristics to simulate the impact of O₃ on crop yield. Model inputs are irradiance, temperature, relative humidity, precipitation, air pressure, wind speed, and O₃ concentration at a reference height (C_z) to calculate atmospheric resistances (R_a) and boundary layer resistances (R_b) for O₃ deposition to the crop canopy. It further incorporates crop-specific parameters related to leaf physiology, phenology, and carbon coefficients, alongside site-specific data (latitude, longitude, and elevation) to simulate crop growth at stages from sowing to maturity, denoted by the development index (DVI). The canopy is divided into four vertical layers, each characterised by the sunlit leaf area index (LAI_sun) and the shaded (LAI_sh) leaf area index, which influence the photosynthetic capacity (V_cmax) and O₃ uptake in each layer. The model accounts for in-canopy resistance (r_inc) and external resistance (r_ext) in each layer, affecting the O₃ flux (accf_st) and its impact on net photosynthesis (A_net) and stomatal conductance ($g_{{\mathrm{O}}_{\mathrm{3}}}$). The A_net–g_sto relationship is modelled using the Leuning (1995) model. Damage from O₃ is estimated after Ewert and Porter (2000) for different canopy layers, which are aggregated to give the overall O₃ impact on canopy A_net, which is integrated according to the JULES-crop model (Osborne et al., 2015), which uses the daily accumulated canopy A_net to calculate the net primary productivity (NPP). The NPP is then distributed as carbon to various parts of the crop (roots (C_root), stems (C_stem), leaves (C_leaf), harvestable organs (C_harv)). The C_harv provides the yield and grain dry matter, C_leaf provides the LAI, and C_stem provides the crop height. The DO₃SE-Crop model requires hourly input meteorological and O₃ concentration data which are used to produce output on either an hourly (i.e. leaf physiology and short-term O₃ damage variables) or daily (i.e. phenology, soil moisture, long-term O₃ damage, C allocation, biomass, and yield variables) time step.

2.2 DO₃SE-Crop phenology

The DO₃SE-Crop model uses thermal time to define the rate of crop development in relation to the timing of three key developmental stages, TT_emr (the period from sowing to emergence), TT_veg (the period of emergence to start of grain filling), and TT_rep (the period from the start of grain filling to maturity), based on the method of Osborne et al. (2015). Thermal time is calculated by accumulating an effective temperature (T_eff) using base (T_b), optimum (T_o), and maximum (T_m) cardinal temperatures as shown in Eq. (1):

$$\begin{array}{}\text{(1)}& T_{\mathrm{eff}}=\left\{\begin{array}{ll}\mathrm{0}& \text{for}{T}_{\mathrm{air}}<T_{\mathrm{b}}\\ T_{\mathrm{air}}-T_{\mathrm{b}}& \text{for}{T}_{\mathrm{b}}\le T_{\mathrm{air}}\le T_{\mathrm{0}}\\ \left(T_{\mathrm{0}}-T_{\mathrm{b}}\right)\left(\mathrm{1}-\frac{T_{\mathrm{air}}-T_{\mathrm{0}}}{T_{\mathrm{m}}-T_{\mathrm{0}}}\right)& \text{for}{T}_{\mathrm{0}}<T_{\mathrm{air}}<T_{\mathrm{m}}\\ \mathrm{0}& \text{for}{T}_{\mathrm{air}}\ge T_{\mathrm{m}}\end{array}\right\},\end{array}$$

where T_air is the surface air temperature in °C, and T_eff is at a maximum when T_air=T_o; this point denotes the highest developmental rate. T_eff declines as the temperature falls or rises above T_o, with a linear decrease in crop development. T_eff is zero, i.e. no development, when T_air falls below or rises above T_b and T_m, respectively; i.e. $T_{\mathrm{m}}\le T_{\mathrm{air}}<T_{\mathrm{b}}$. During the sowing to emergence phase, development is dependent on T_b, whereas during the vegetative and reproductive phase, development depends on T_m or T_o.

Winter wheat requires vernalisation (a period of exposure to low temperature during germination to accelerate flowering). Vernalisation alters the length of TT_veg and hence flowering initiation, with subsequent effects on later growth stages such as heading. Vernalisation occurs when the minimum (VT_min) and maximum (VT_max) daily temperature is less than 15 and 30 °C, respectively (Zheng et al., 2015). Accumulated vernalised days (V_dd) are calculated as the sum of vernalised and devernalised days from emergence to the start of anthesis (Zheng et al., 2015) as shown in Eq. (2):

$$\begin{array}{}\text{(2)}& V_{dd}=\sum (V-V_d),\end{array}$$

where $V=\left(\mathrm{1.4}-\mathrm{0.778}\times T_{\mathrm{air}},\mathrm{0.5}+\mathrm{13.44}\frac{T_{\mathrm{air}}}{{\left(T_{\max }-T_{\min }+\mathrm{3}\right)}^{\mathrm{2}}}\right)$ for VT_max <30 °C and VT_min<15 °C, $V_d=(min\left(\mathrm{0.5}\left(T_{\max }-\mathrm{30}\right),V_{\mathrm{prev}}\right)$ for VT_max>30 °C, and V_dd<10 d.

The vernalisation factor (VF) decreases from 1 to 0 as (V_dd) increases. VF depends on a cultivar-specific vernalisation coefficient (PIV) as described by Eq. (3):

$$\begin{array}{}\text{(3)}& \mathrm{VF}=\mathrm{1}-(\mathrm{0.0054545}\times \mathrm{PIV}+\mathrm{0.0003})\times (\mathrm{50}-V_{dd}).\end{array}$$

Photoperiod (PP) or day length also affects the occurrence and timing of the flowering stage and is calculated according to latitude using standard solar geometry to estimate day length (Jones, 1992). The photoperiod factor (PF) represents the sensitivity to PP, which decreases from 1 to 0 as the photoperiod shortens and is estimated according to a cultivar-specific photoperiod coefficient (PID) after Tao et al. (2012) as described in Eq. (4):

$$\begin{array}{}\text{(4)}& \mathrm{PF}=\mathrm{1}-\left[\left({\displaystyle \frac{\mathrm{PID}}{\mathrm{10}\mathrm{000}}}\right)\times {\left(\mathrm{20}-\mathrm{PP}\right)}^{\mathrm{2}}\right].\end{array}$$

Crop development is related to the development index (DVI) after Osborne et al. (2015), which takes values of −1 upon sowing, 0 on emergence, 1 at anthesis, and 2 at crop maturity. The DO₃SE-Crop model DVI equations have been modified from Osborne et al. (2015) to take account of the photoperiod and vernalisation for winter wheat (see Eq. 5); for spring wheat these factors are omitted:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& -\mathrm{1}\le \mathrm{DVI}<\mathrm{0}\mathrm{for}{\mathrm{TT}}_{\mathrm{eff}}<{\mathrm{TT}}_{\mathrm{emr}}\\ & \mathrm{0}\le \mathrm{DVI}<\mathrm{1}\mathrm{for}{\mathrm{TT}}_{\mathrm{emr}}\le {\mathrm{TT}}_{\mathrm{eff}}\times \mathrm{VF}\times \mathrm{PF}<{\mathrm{TT}}_{\mathrm{veg}}\\ & \mathrm{1}\le \mathrm{DVI}\le \mathrm{2}\mathrm{for}{\mathrm{TT}}_{\mathrm{veg}}\le {\mathrm{TT}}_{\mathrm{eff}}\le {\mathrm{TT}}_{\mathrm{rep}}.\end{array}\end{array}$$

DO₃SE-Crop allows for any number of representative leaf populations (pop) and canopy layers (n) to be defined over the course of the crop growing season by dividing leaf populations as they emerge evenly across the canopy layers defined by LAI. In this study, we used a single leaf population and four canopy layers (i.e. pop =1; n=4) for simplicity. The crop sowing is assumed to be at DVI $=-\mathrm{1}$ (start of TT_emr) and emergence at DVI =0 (start of TT_veg). The flag leaf is assumed to develop at DVI =1, at the commencement of TT_rep, marking the initiation of anthesis (A_start flowering) and flag leaf emergence, which typically occurs 4–5 d prior to the onset of anthesis and is further divided into expanding and senescing leaf periods (i.e. tl_ep and tl_se) with a default ratio of 0.67 to 0.33 for each of these periods. Maturity is assumed at DVI =2, at the end of TT_rep. The model allows estimation of the POD_y metric by accumulating stomatal O₃ flux from the start of anthesis to maturity. The total canopy leaf life span (TT_leaf) of the crop is distributed over the DVI between 0 and 2. The total lifespan (T_l) covers the full period from sowing to maturity, corresponding to DVI between −1 and 2. The relationship between these different variables is described in Fig. 2.

2.3 DO₃SE-Crop leaf-level physiology

Key leaf-level physiological variables of the DO₃SE-Crop model are A_net and g_sto. Net photosynthesis is simulated using the biochemical photosynthesis-based model initially developed by Farquhar et al. (1980) and since modified by Sharkey et al. (2007). The coupled A_net–g_sto model of Leuning (1995) is used to estimate g_sto from A_net, which means that g_sto is regulated by the demand of CO₂ for A_net on consideration of environmental conditions and crop physiology. Ozone stress, causing both instantaneous effects on A_net and long-term effects on A_net via leaf senescence, is simulated based on algorithms developed by Ewert and Porter (2000).

2.3.1 Leaf net photosynthesis (A_net)

The A_net model assumes that photosynthesis is constrained depending on prevailing environmental conditions according to three main mechanisms: rubisco activity (A_c); ribulose-1,5-bisphosphate (RuBP) regeneration, which is constrained by the speed of electron transport (A_j); and the low rate of transfer of photosynthetic products (most frequently triose phosphate consumption) (A_p) (Sharkey et al., 2007) and by soil water stress (f_PAW). The algorithm for A_c, which is based on Medlyn et al. (2002) and modified in DO₃SE-Crop to include the O₃ damage functions, is given in Eq. (6):

$$\begin{array}{}\text{(6)}& A_c=V_{\mathrm{cmax}}\times f_{\mathrm{PAW}}\times {\displaystyle \frac{\left(C_{\mathrm{i}}-{\mathrm{\Gamma }}^{\ast }\right)\times f{\mathrm{O}}_{\mathrm{3},\mathrm{s}}\left(d\right)\times f_{\mathrm{LS}}}{C_{\mathrm{i}}+K_{\mathrm{c}}\left(\mathrm{1}+\frac{O_{\mathrm{i}}}{K_o}\right)}},\end{array}$$

where V_cmax (µmol CO₂ m⁻² s⁻¹) is the maximum carboxylation capacity at 25 °C, C_i (µmol mol⁻¹) and O_i (mmol mol⁻¹) are the intercellular CO₂ and O₂ partial pressures, K_c (µmol mol⁻¹) and K₀ (mmol mol⁻¹) are the rubisco Michaelis–Menten constants for CO₂ and O₂, Γ^∗ (µmol mol⁻¹) is the CO₂ compensation point in the absence of respiration, fO_3,s(d) is the factor that accounts for the cumulative stomatal O₃ flux effect on V_cmax over the course of a day, and f_LS is the factor that accounts for the cumulative stomatal O₃ flux effect over the course of a leaf life span on leaf senescence. Section 2.3.2 gives a full description of the methods used to estimate O₃ damage. The f_PAW factor is calculated by Eq. (7):

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& f_{\mathrm{PAW}}=\mathrm{1}{\text{for\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}PAW}}_t\le \mathrm{PAW}\le \mathrm{100}\mathit{\%},\\ & f_{\mathrm{PAW}}=\mathrm{1}+\mathit{\left\{}{\displaystyle \frac{\mathrm{PAW}-{\mathrm{PAW}}_t}{{\mathrm{PAW}}_t}}\mathit{\right\}}\text{for\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}PAW}\le {\mathrm{PAW}}_t.\end{array}\end{array}$$

PAW is the amount of water in the soil (in % terms) which is available to the plant estimated according to the DO₃SE models' single-soil-layer bucket model (Büker et al., 2012). At PA =100 % the soil is at field capacity; at PAW $=\mathrm{0},\mathit{\%}$ the soil is at wilting point. PAW_t is the threshold PAW, above which it is assumed there is no constraint on A_c, defined as 50 % after LRTAP (2017). Only once PAW<PAW_t will soil water begin to limit g_sto and hence stomatal O₃ flux.

The constraint on photosynthesis due to the rate of electron transport A_j is described in Eq. (8):

$$\begin{array}{}\text{(8)}& A_j=J\times {\displaystyle \frac{C_{\mathrm{i}}-{\mathrm{\Gamma }}^{\ast }}{a\times C_{\mathrm{i}}+b\times {\mathrm{\Gamma }}^{\ast }}},\end{array}$$

where J is the electron transport rate (µmol CO₂ m⁻² s⁻¹), and the parameters a and b denote the electron requirements for the formation of NADPH and ATP, respectively (Sharkey et al., 2007)

Finally, the photosynthesis limitation due to the low rate of transfer of photosynthetic products A_p (µmol CO₂ m⁻² s⁻¹) is given in Eq. (9):

$$\begin{array}{}\text{(9)}& A_{\mathrm{p}}=\mathrm{0.5}\times V_{\mathrm{cmax}}.\end{array}$$

The leaf net photosynthesis (A_net) in µmol CO₂ m⁻² s⁻¹ is calculated by Eq. (10):

$$\begin{array}{}\text{(10)}& A_{\mathrm{net}}=\left(A_c,A_j,A_{\mathrm{p}}\right)-R_{\mathrm{d}},\end{array}$$

where leaf dark respiration (R_d) in µmol CO₂ m⁻² s⁻¹ is calculated as Vc_max×R_dcoeff, where R_dcoeff is the leaf dark respiration coefficient initially set equal to 0.015 after Clark et al. (2011) – a value provided for C₃ grasses.

2.3.2 Short- and long-term O₃ damage to A_c

The short-term impact of O₃ on A_c is calculated according to the fO_3,s(d) factor (between 0 and 1), which allows for an instantaneous effect of O₃ on photosynthesis when stomatal O₃ flux (f_st), in nmol O₃ m⁻² s⁻¹ and calculated as described later in Sect. 1.2.3, overwhelms detoxification and repair mechanisms (Betzelberger et al., 2012; Y. Feng et al., 2022) and is estimated following Ewert and Porter (2000). Here, $f_{{\mathrm{O}}_{\mathrm{3}},\mathrm{s}}\left(h\right)$ represents the relationship between f_st and a potential decrease in A_c calculated for every hour of the day by Eq. (11):

$$\begin{array}{}\text{(11)}& \begin{array}{rl}& f_{{\mathrm{O}}_{\mathrm{3}},\mathrm{s}}\left(h\right)=\mathrm{1};\text{for}{f}_{\mathrm{st}}\le {\displaystyle \frac{\mathit{\gamma }\mathrm{1}}{\mathit{\gamma }\mathrm{2}}}\\ & f_{{\mathrm{O}}_{\mathrm{3}},\mathrm{s}}\left(h\right)=\mathrm{1}+\mathit{\gamma }\mathrm{1}-\mathit{\gamma }\mathrm{2}\times f_{\mathrm{st}}\\ & \mathrm{for}{\displaystyle \frac{\mathit{\gamma }\mathrm{1}}{\mathit{\gamma }\mathrm{2}}}<f_{\mathrm{st}}<{\displaystyle \frac{\mathrm{1}+\mathit{\gamma }\mathrm{1}}{\mathit{\gamma }\mathrm{2}}}\\ & f_{{\mathrm{O}}_{\mathrm{3}},\mathrm{s}}\left(h\right)=\mathrm{0};\mathrm{for}{f}_{\mathrm{st}}\ge {\displaystyle \frac{\mathrm{1}+\mathit{\gamma }\mathrm{1}}{\mathit{\gamma }\mathrm{2}}},\end{array}\end{array}$$

where γ1 (dimensionless) and γ2 (nmol O₃ m⁻² s⁻¹)⁻¹ are both short-term O₃ damage coefficients, with γ1 representing the O₃ detoxification threshold below which no damage occurs to the photosynthetic system and γ2 determines the effect of f_st on A_c once this detoxification threshold is exceeded; fO_3,s(d) and $f{\mathrm{O}}_{\mathrm{3},\mathrm{s}}\left(d-\mathrm{1}\right)$ (i.e. fO_3,s(d) at the end of the previous day) are calculated by Eq. (12),

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{f}_{{\mathrm{O}}_{\mathrm{3}},\mathrm{s}}\left(d\right)& =f_{{\mathrm{O}}_{\mathrm{3}},\mathrm{s}}\left(h\right)\times r_{{\mathrm{O}}_{\mathrm{3}},\mathrm{s}};\\ & \text{for}\mathrm{PAR}\le \mathrm{50}{\mathrm{W}\mathrm{m}}^{-\mathrm{2}}\\ f_{{\mathrm{O}}_{\mathrm{3}},\mathrm{s}}\left(d\right)& =f_{{\mathrm{O}}_{\mathrm{3}},\mathrm{s}}\left(h\right)\times f_{{\mathrm{O}}_{\mathrm{3}},\mathrm{s}}(d-\mathrm{1})\\ & \text{for}\mathrm{PAR}>\mathrm{50}{\mathrm{W}\mathrm{m}}^{-\mathrm{2}},\end{array}\end{array}$$

where $r_{{\mathrm{O}}_{\mathrm{3}},\mathrm{s}}$ (dimensionless) represents incomplete recovery from O₃ overnight, which depends on leaf age according to Eq. (13),

$$\begin{array}{}\text{(13)}& r_{{\mathrm{O}}_{\mathrm{3}},\mathrm{s}}=f_{{\mathrm{O}}_{\mathrm{3}},\mathrm{s}}(d-\mathrm{1})+\left(\mathrm{1}-f_{{\mathrm{O}}_{\mathrm{3}},\mathrm{s}}(d-\mathrm{1})\right)\times f_{\mathrm{LA}}.\end{array}$$

The long-term impact of O₃ on V_cmax represented by the f_Ls term represents the longer-term accumulation of stomatal O₃ flux (acc_fst), causing degradation to the rubisco enzyme, which triggers early and enhanced senescence of mature leaves (Gelang et al., 2000; Osborne et al., 2019). The acc_fst term is accumulated from 200 °C days before anthesis until maturity to be consistent with LRTAP (2017), which defines this as the O₃-sensitive period for wheat. The simulation of f_Ls (and f_LA used in the short-term O₃ effect) is related to thermal-time-defined periods over the course of a leaf population life span TT_leaf as described in Fig. 2.

Figure 2 The division of thermal-time-defined periods (TT_emr, TT_veg, TT_rep, and TT_leaf and the relationship with f_LA and f_Ls) for the canopy, as represented in this study by a single leaf population.

The O₃ effect on f_Ls is first simulated by estimating a weighted accumulated fst ($f{{\mathrm{O}}_{\mathrm{3}}}_l$) modified from Ewert and Porter (2000) by Eq. (14):

$$\begin{array}{}\text{(14)}& f{{\mathrm{O}}_{\mathrm{3}}}_l=\mathrm{1}-\max (\min \left(\mathit{\gamma }\mathrm{3}\times \left({\mathrm{acc}}_{\mathrm{fst}}-\mathrm{CLsO}\mathrm{3}\right),\mathrm{1}\right),\mathrm{0}),\end{array}$$

where γ3 determines the occurrence of senescence once a critical cumulative stomatal O₃ flux CLsO3 (in mmol m⁻²) has been exceeded. The rate of senescence is determined by γ4, which determines the onset of senescence, and γ5, which determines maturity as described in Eq. (15):

$$\begin{array}{}\text{(15)}& \begin{array}{rl}& t{l}_{{\mathrm{ep}}_{{\mathrm{O}}_{\mathrm{3}}}}=t{l}_{\mathrm{ep}}\times (\mathrm{1}-\left(\left(\mathrm{1}-f{{\mathrm{O}}_{\mathrm{3}}}_l\right)\times \mathit{\gamma }\mathrm{4}\right))\\ & t{l}_{{\mathrm{se}}_{{\mathrm{O}}_{\mathrm{3}}}}=t{l}_{\mathrm{se}}\times \left(\mathrm{1}-\left(\left(\mathrm{1}-f{{\mathrm{O}}_{\mathrm{3}}}_l\right)\times \mathit{\gamma }\mathrm{5}\right)\right)+zc\\ & zc=t{l}_{\mathrm{ep}}-t{l}_{{\mathrm{ep}}_{{\mathrm{O}}_{\mathrm{3}}}},\end{array}\end{array}$$

where tl_ep is the thermal time accumulated by a leaf (LTT) in °C days between a fully expanded leaf and the start of leaf senescence, $t{l}_{{\mathrm{ep}}_{{\mathrm{O}}_{\mathrm{3}}}}$ is tl_ep with an O₃ effect which may bring senescence earlier, tl_se is the LTT in °C days between the onset of senescence and maturity, and $t{l}_{{\mathrm{se}}_{{\mathrm{O}}_{\mathrm{3}}}}$ is tl_se with an O₃ effect which may bring maturity earlier. f_Ls is estimated by Eq. (16):

$$\begin{array}{}\text{(16)}& \begin{array}{rl}{f}_{\mathrm{Ls}}& =\mathrm{1};\text{for}\mathrm{LTT}\le T{T}_{\mathrm{veg}}+t{l}_{\mathrm{ep}}\\ f_{\mathrm{Ls}}& =\mathrm{1}-{\displaystyle \frac{\mathrm{LTT}-T{T}_{\mathrm{veg}}-t{l}_{{\mathrm{ep}}_{{\mathrm{O}}_{\mathrm{3}}}}}{t{l}_{{\mathrm{se}}_{{\mathrm{O}}_{\mathrm{3}}}}}};\\ & \text{for}T{T}_{\mathrm{veg}}+t{l}_{\mathrm{ep}}<\mathrm{LTT}<T{T}_{\mathrm{leaf}}\\ f_{\mathrm{Ls}}& =\mathrm{0};\text{for}\mathrm{LTT}\ge T{T}_{\mathrm{leaf}}.\end{array}\end{array}$$

2.3.3 Stomatal conductance (g_sto)

The coupled photosynthesis–stomatal-conductance (A_net–g_sto) model based on Leuning (1995) and modified for vapour pressure deficit (VPD) is used to estimate $g_{{\mathrm{CO}}_{\mathrm{2}}}$ and stomatal conductance to CO₂ in µmol CO₂ m⁻² s⁻¹ as described in Eq. (17),

$$\begin{array}{}\text{(17)}& g_{{\mathrm{CO}}_{\mathrm{2}}}=[f_{\min }+m\times A_{\mathrm{net}}\times f_{\mathrm{VPD}}/(c_{\mathrm{s}}-\mathrm{\Gamma }\left)\right],\end{array}$$

where f_min (µmol m⁻² s⁻¹) is the minimum daytime $g_{{\mathrm{CO}}_{\mathrm{2}}}$ (Leuning, 1990). The parameter m (dimensionless) is the composite sensitivity of $g_{{\mathrm{CO}}_{\mathrm{2}}}$ to assimilation rate and vapour pressure deficit (VPD), with the relationship between VPD and relative stomatal conductance (f_VPD) estimated by Eq. (18),

$$\begin{array}{}\text{(18)}& f_{\mathrm{VPD}}={\left(\mathrm{1}+{\left({\displaystyle \frac{\mathrm{VPD}}{\mathrm{VPDo}}}\right)}^{\mathrm{8}}\right)}^{-\mathrm{1}},\end{array}$$

where VPD₀ is an empirical parameter, defined using boundary line analysis, describing the variation in relative stomatal conductance with VPD (Danielsson et al., 2003; Pleijel et al., 2007). c_s (mmol mol⁻¹) is the external CO₂ concentration at the leaf surface and is calculated from the external CO₂ concentration at the upper surface of the leaf boundary layer c_a (mmol mol⁻¹) so that $c_{\mathrm{s}}=c_a-\left(\frac{A_{\mathrm{net}}}{g_{{\mathrm{bCO}}_{\mathrm{2}}}}\right)$ after Masutomi (2023), where $g_{{\mathrm{bCO}}_{\mathrm{2}}}$ is the boundary layer conductance to CO₂ (in mol m⁻² s⁻¹), and conversion factors for gases and heat across the boundary layer are given in Sect. S1a.

Finally, $g_{{\mathrm{CO}}_{\mathrm{2}}}$ is converted to $g_{{\mathrm{O}}_{\mathrm{3}}}$ in mmol O₃ m⁻² s⁻¹ by dividing by 1000 and using the conversion factor 0.96, which assumes that the ratio of the diffusivities of gases in air is equal to the inverse of the square root of the ratio of molecular weights (as described in Campbell and Norman, 1998); see also Sect. S1b).

2.3.4 Stomatal ozone flux (f_st)

Stomatal [O₃] flux (f_st in nmol m⁻² s⁻¹) is calculated after the method described in the UNECE mapping manual (LRTAP, 2017) described in Eq. (19):

$$\begin{array}{}\text{(19)}& f_{\mathrm{st}}=C_{\mathrm{l}}\times g_{{\mathrm{O}}_{\mathrm{3}}\mathrm{m}/\mathrm{s}}\times {\displaystyle \frac{r_{\mathrm{c}}}{r_{\mathrm{b},{\mathrm{O}}_{\mathrm{3}}}+r_{\mathrm{c}}}},\end{array}$$

where C_l is the [O₃] at the upper surface of the laminar layer of a leaf (nmol O₃ m⁻³). Ozone concentration in ppb can be converted to nmol m⁻³ by multiplying O₃ in ppb by P/($R\times T_{\mathrm{air},\mathrm{k}})$, where P is the atmospheric pressure (1.013×10⁵ in Pascal), R is the universal gas constant (8.31447 J mol K⁻¹), and T_air,k is surface air temperature in °K. To convert $g_{{\mathrm{O}}_{\mathrm{3}}}$ (mol O₃ m⁻² s⁻¹) to $g_{{\mathrm{O}}_{\mathrm{3}}\mathrm{m}/\mathrm{s}}$ (m s⁻¹) we assume a standard temperature (20 °C) and P divided by 41 to give the conductance value in m s⁻¹. The $r_{\mathrm{c}}/(r_{\mathrm{b},{\mathrm{O}}_{\mathrm{3}}}+r_{\mathrm{c}})$ term represents the O₃ deposition rate to the leaf through resistances r_b (the quasi-laminar resistance (s m⁻¹)) and r_c (the leaf surface resistance (s m⁻¹)), which allow for both stomatal and non-stomatal deposition to the leaf surface. r_c is $\mathrm{1}/(g_{{\mathrm{O}}_{\mathrm{3}}\mathrm{m}/\mathrm{s}}+g_{\mathrm{ext}})$, where g_ext is $\mathrm{1}/\mathrm{2500}$ (s m⁻¹). $r_{\mathrm{b},{\mathrm{O}}_{\mathrm{3}}}$ is estimated by Eq. (20):

$$\begin{array}{}\text{(20)}& r_{\mathrm{b},{\mathrm{O}}_{\mathrm{3}}}=\mathrm{1.3}\times \mathrm{150}\times \sqrt{{\displaystyle \frac{L}{u_{\mathrm{l}}}}},\end{array}$$

where the factor 1.3 accounts for the differences in diffusivity between heat and O₃ (see Sect. S1a), L is the cross-wind leaf dimension (m), and u_l is the wind speed (m s⁻¹) at the top of the leaf laminar boundary layer. The leaf boundary layer resistance to CO₂ is estimated using a value of 1.24 for the difference between heat and CO₂ in place of the 1.3 value for O₃ (Campbell and Norman, 1998).

2.4 DO₃SE-Crop canopy

The DO₃SE-Crop model uses a multi-layer approach to scale from leaf to the canopy. We assume that wind, irradiance, [O₃] concentration, and leaf nitrogen content are the key environmental conditions which change with the cumulative canopy leaf area index (LAI) and influence leaf physiology and therefore canopy layer estimates of A_net, $g_{{\mathrm{O}}_{\mathrm{3}}}$, and g_ext; other environmental variables (e.g. T_air and VPD) are assumed to remain constant over the canopy.

2.4.1 Canopy irradiance

Changes in irradiance through the canopy are described as sunlit and shaded canopy fractions and the associated quantity of direct and diffuse photosynthetically active radiation (PAR, W m⁻²); these are estimated according to increasing levels of cumulative LAI using the methods of (de Pury and Farquhar, 1997); full details are given in the Sect. S2. Application of this method requires the canopy to be divided into layers of equal LAI, including both green LAI (LAI_G) and brown (LAI_B) LAI.

PAR absorbed per unit leaf area is divided into PAR_dir and PAR_diff, which also includes scattered (re-reflected by the canopy) beam calculated by

$$\begin{array}{}\text{(21)}& {\displaystyle }{\mathrm{PAR}}_{\mathrm{dir}}\left(\mathrm{LAI}\right)=(\mathrm{1}-{\mathit{\rho }}_{\mathrm{cb}}(\mathit{\beta }\left)\right)k_{\mathrm{b}}{}^{\prime }{I}_{\mathrm{b}}\left(\mathrm{0}\right)\exp (-k_{\mathrm{b}}^{\prime }\mathrm{LAI}),\text{(22)}& {\displaystyle }{\mathrm{PAR}}_{\mathrm{diff}}\left(\mathrm{LAI}\right)=\left(\mathrm{1}-{\mathit{\rho }}_{\mathrm{cd}}\right)k_{\mathrm{d}}{I}_{\mathrm{d}}\left(\mathrm{0}\right)\exp (-k_{\mathrm{d}}\mathrm{LAI}),\end{array}$$

where PAR_dir is the absorbed beam plus scattered beam PAR (photosynthetically active radiation) per unit leaf area, PAR_diff is the absorbed diffuse plus scattered diffuse PAR per unit leaf area, ρ_cb is the canopy reflection coefficient for beam PAR, ρ_cd is the canopy reflection coefficient for diffuse PAR, $k_{\mathrm{b}}^{\prime }$ is the beam and scattered beam PAR extinction coefficient, $k_{\mathrm{d}}^{\prime }$ is the diffuse and scattered diffuse PAR extinction coefficient, β is the solar elevation above the horizontal plane of the Earth's surface, I_b(0) is beam PAR per unit ground area at the top of the canopy, and I_d(0) is diffuse PAR per unit ground area at the top of the canopy.

Estimates of the LAI fractions of sunlit (LAI_sun) and shaded (LAI_sh) parts of each canopy layer (i) are made by Eqs. (23) and (24):

$$\begin{array}{}\text{(23)}& {\mathrm{LAI}}_{\mathrm{sun},i}=\left[\mathrm{1}-\exp \left(-\mathrm{0.5}\times {\displaystyle \frac{\mathrm{LAI}i}{\sin \mathit{\beta }}}\right)\right]\times \mathrm{2}\sin \mathit{\beta },\end{array}$$

where β is the solar elevation angle (see Sect. S3), and

$$\begin{array}{}\text{(24)}& {\mathrm{LAI}}_{\mathrm{sh},i}={\mathrm{LAI}}_i-{\mathrm{LAI}}_{\mathrm{sun},i}.\end{array}$$

The DO₃SE-Crop model simulates LAI as part of the crop growth model, and LAI is assumed to be evenly distributed across all layers (see Sect. 1.4.2 and Eq. 43).

Therefore, PAR for the sunlit part of each layer (PAR_sun) can be described as

$$\begin{array}{rl}{\int }_{{\mathrm{LAI}}_i}^{{\mathrm{LAI}}_n}{\mathrm{PAR}}_{\mathrm{sun}}& ={\int }_{{\mathrm{LAI}}_i}^{{\mathrm{LAI}}_n}\left({\mathrm{LAI}}_{\mathrm{sun},i}\right)\times ({\mathrm{PAR}}_{\mathrm{sh}}\\ & +{\mathrm{PAR}}_{\mathrm{bsun}}\left(\mathit{\beta }\right))\mathrm{dLAI},\end{array}$$

where PAR_sh is absorbed PAR by shaded leaves per unit leaf area; PAR_bsun is beam PAR absorbed by sunlit leaves per unit leaf area; and ${\int }_{{\mathrm{LAI}}_i}^{{\mathrm{LAi}}_n}{\mathrm{PAR}}_{\mathrm{dir}}$ can be written as $(\mathrm{1}-\mathit{\rho }cb(\mathit{\beta }\left)\right)\times k_{\mathrm{b}}{}^{\prime }\times Ib\times \left[\exp \right(-k_{\mathrm{b}}{}^{\prime }{\mathrm{LAI}}_i)-\exp (-k_{\mathrm{b}}{}^{\prime }{\mathrm{LAI}}_n\left)\right]$ and ${\mathrm{PAR}}_{\mathrm{bsun}}\left(\mathit{\beta }\right)=\left(\mathrm{1}-\mathit{\sigma }\right)I_{\mathrm{b}}\left(\mathrm{0}\right)\frac{{\mathrm{Cos}}_{{\mathit{\alpha }}_{\mathrm{l}}}}{\mathrm{Sin}\mathit{\beta }}$, where α_l is angle of irradiance beam on the leaf normal and σ is the leaf scattering coefficient for PAR.

Similarly, PAR for the shaded part of each layer (PAR_sun) can be described as

$${\int }_{{\mathrm{LAI}}_i}^{{\mathrm{LAI}}_n}{\mathrm{PAR}}_{\mathrm{sh}}={\int }_{{\mathrm{LAI}}_i}^{{\mathrm{LAI}}_n}\left({\mathrm{LAI}}_{\mathrm{sh},i}\right)\times ({\mathrm{PAR}}_{\mathrm{diff}}+{\mathrm{PAR}}_{\mathrm{bsun}})\mathrm{dLAI},$$

where ${\int }_{{\mathrm{LAi}}_i}^{{\mathrm{LAi}}_n}\left({\mathrm{PAR}}_{\mathrm{diff}}\right(\mathrm{LAI})$ can be written as $(\mathrm{1}-{\mathit{\rho }}_{\mathrm{cd}})\times k_{\mathrm{b}}{}^{\prime }\times Ib\times \left[\exp \right(-k_{\mathrm{d}}^{\prime }{\mathrm{LAI}}_i)-\exp (-k_{\mathrm{d}}^{\prime }{\mathrm{LAI}}_n\left)\right]\mathrm{d}L$, ${\int }_{{\mathrm{LAI}}_i}^{{\mathrm{LAi}}_n}{\mathrm{PAR}}_{\mathrm{bs}}\left(\mathrm{LAI}\right)$ is $I_{\mathrm{b}}\left(\mathrm{0}\right)[{\mathrm{PAR}}_{\mathrm{dir}}-\left(\mathrm{1}-\mathit{\sigma }\right)k_{\mathrm{b}}\times [\exp \left(-k_{\mathrm{b}}{\mathrm{LAI}}_i\right)-\exp \left(-k_{\mathrm{b}}{\mathrm{LAI}}_n\right)]$, and PAR_bs (LAI) is absorbed scattered beam PAR per unit leaf area.

2.4.2 Canopy [O₃] concentration

O₃ concentration will vary as a function of O₃ loss to the canopy (i.e. deposition via the stomates and external plant parts) and O₃ replacement from ambient air concentrations above the canopy. Limited data have been collected showing how O₃ concentrations vary with canopy depth in semi-natural communities (Jaggi et al., 2006). These data suggest that a minimum, bottom of canopy O₃ concentration (C_zb) is about 0.2 times that at the top of the canopy (C_zh) and that the O₃ concentration difference within the canopy is closely related to the LAI of the canopy layers.

Since each canopy layer can be assumed to be a parallel sink, the O₃ flux to a layer depends on the conductance (inverse of resistance) of that layer and the O₃ concentration at the top of the layer (C_i; with C₀ being C_zh (i.e. the O₃ concentration at height C_h, the top of the canopy)); we follow and generalise the work of Waggoner (1971) by separating the canopy into nL leaf layers. We calculate the O₃ concentration for each layer, C_i, from O₃ intake, I_i, by

$$\begin{array}{}\text{(25)}& C_{\mathrm{i}}=r_{c,i}{I}_i,\end{array}$$

with r_c,i the leaf surface resistance to O₃ for layer i. I_i is calculated as the solution to a system of linear equations. The in-canopy aerodynamic resistance for layer i is described in terms of r_c,i, I_i, and the resistances of the bulk air among the leaves (R_i). Assuming a uniform O₃ concentration C₀ above the canopy, we use generalised equations from Waggoner (1971) to calculate the difference in O₃ concentration between the exterior air and the leaf interior. For the top layer, this difference is C₀ minus 0, while for each lower layer, the difference decreases progressively depending on the resistances and fluxes within the canopy. This O₃ concentration difference is calculated by

$$\begin{array}{}\text{(26)}& C_{\mathrm{0}}=R_i\sum _{j=\mathrm{1}}^{\mathrm{nL}}{I}_j+r_{c,\mathrm{1}}{I}_{\mathrm{1}}\end{array}$$

for the top canopy layer,

$$\begin{array}{}\text{(27)}& \mathrm{0}=R_i\sum _{j=i}^{\mathrm{nL}}{I}_j+r_{c,i}{I}_i-r_{c,i-\mathrm{1}}{I}_{i-\mathrm{1}}\end{array}$$

for each canopy layer i between the top layer and the bottom layer, and

$$\begin{array}{}\text{(28)}& \mathrm{0}=R_{\mathrm{nL}+\mathrm{1}}{I}_{\mathrm{nL}+\mathrm{1}}-r_{c,\mathrm{nL}}{I}_{\mathrm{nL}}\end{array}$$

for the bottom layer of the canopy between the lowest leaf layer and the ground. These can also be written into the matrix form

$$\begin{array}{}\text{(29)}& \begin{array}{rl}& \left(\begin{array}{ccccc}{r}_{c,\mathrm{1}}+R_{\mathrm{1}}& R_{\mathrm{1}}& R_{\mathrm{1}}& \mathrm{\cdots }& R_{\mathrm{1}}\\ -r_{c,\mathrm{1}}& r_{c,\mathrm{2}}+R_{\mathrm{2}}& R_{\mathrm{2}}& \mathrm{\cdots }& R_{\mathrm{2}}\\ \mathrm{0}& -r_{c,\mathrm{2}}& r_{c,\mathrm{3}}+R_{\mathrm{3}}& \mathrm{\cdots }& R_{\mathrm{3}}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{⋮}& R_{nL+\mathit{\_}\mathrm{1}}\end{array}\right)\\ & \left(\begin{array}{c}{I}_{\mathrm{1}}\\ I_{\mathrm{2}}\\ I_{\mathrm{3}}\\ \mathrm{⋮}\\ I_{nL+\mathrm{1}}\end{array}\right)=\left(\begin{array}{c}{C}_{\mathrm{0}}\\ \mathrm{0}\\ \mathrm{0}\\ \mathrm{⋮}\\ \mathrm{0}\end{array}\right),\end{array}\end{array}$$

which can be numerically solved for I_x when $r_{c,\mathrm{1}}\ne \mathrm{0}$ and R₁≠0.

Resistances for each layer are calculated as described in the Supplement (Sect. S5) using standard DO₃SE deposition modelling methods (Emberson et al., 2000, 2001).

2.4.3 Canopy maximum carboxylation capacity (V_cmax)

We allow for an exponential decrease in leaf N with canopy depth, which will influence both the photosynthetic capacity (V_cmax) and dark respiration (R_dc). Photosynthetic capacity at each canopy layer i is calculated by Eq. (30):

$$\begin{array}{}\text{(30)}& V_{\mathrm{cmax},i}=n_e\times n_{\mathrm{0}}\times e^{-kN\left(\frac{\mathrm{LAI}i}{\mathrm{LAI}}\right)},\end{array}$$

where n_e (mol CO₂ m⁻² s⁻¹ kg C (kg N)⁻¹) is a constant relating leaf nitrogen to rubisco carboxylation capacity, n₀ (kg N[kg C]⁻¹) is the leaf N concentration at the top of the canopy, and kN is a nitrogen profile coefficient initially set at 0.78 after (Clark et al., 2011). The model assumes non-limiting conditions for soil nitrogen, in accordance with the experimental data.

2.4.4 Canopy photosynthesis (Anet_c)

Net canopy photosynthesis (Anet_c) determines the amount of C assimilated by the entire canopy that can subsequently be allocated to different plant parts (i.e. less than the C respired for plant growth and maintenance; see Sect. 1.4.1); the amount of C assimilation will ultimately determine whole-plant biomass. The net photosynthesis for each canopy layer (Anet_i) is calculated according to the LAI fraction of that layer that is sunlit (LAI_sun,i) and shaded (LAI_sh,i) within the layer (i), multiplied by the net photosynthesis of the sunlit (Anet${}_{\mathrm{sun},i,j}$) and shaded leaf (Anet${}_{\mathrm{sh},i,j}$), respectively, described by Eqs. (31) and (32):

$$\begin{array}{}\text{(31)}& {\mathrm{Anet}}_i={\mathrm{LAI}}_{\mathrm{sun},i}\times {\mathrm{Anet}}_{\mathrm{sun},i}+{\mathrm{LAI}}_{\mathrm{sh},i}\times {\mathrm{Anet}}_{\mathrm{sh},i},\end{array}$$

with Anet_c calculated by

$$\begin{array}{}\text{(32)}& {\mathrm{Anet}}_{\mathrm{c}}={\sum }_{i=\mathrm{1}}^n{\mathrm{Anet}}_i.\end{array}$$

Anet_c is converted from µmol CO₂ m⁻² s⁻¹ to kg C m⁻² d⁻¹ by multiplying by 3600 (converting from seconds to hours), multiplying by 1.2 (representing the kg of C per mol), and summing each hourly Anet_c over the course of a day. This Anet_c is used in the Eq. (37).

2.4.5 Canopy stomatal conductance ($g_{{{\mathrm{O}}_{\mathrm{3}}}_c}$)

Similarly, canopy layer (i) stomatal conductance to O₃ ($g_{{{\mathrm{O}}_{\mathrm{3}}}_i}$) is converted from $g_{{\mathrm{CO}}_{\mathrm{2}}}$ by assuming a diffusivity ratio of 0.96 to convert from CO₂ to O₃ and is calculated by Eq. (33) with whole-canopy stomatal conductance calculated by Eq. (34):

$$\begin{array}{}\text{(33)}& {\displaystyle }{g}_{{{\mathrm{O}}_{\mathrm{3}}}_i}={\mathrm{LAI}}_{\mathrm{sun},i}\times g_{{{\mathrm{O}}_{\mathrm{3}}}_{\mathrm{sun},i}}+{\mathrm{LAI}}_{\mathrm{sh},i}\times g_{{{\mathrm{O}}_{\mathrm{3}}}_{\mathrm{sh},i}},\text{(34)}& {\displaystyle }{g}_{{{\mathrm{O}}_{\mathrm{3}}}_c}={\sum }_{i=\mathrm{1}}^n{g}_{{{\mathrm{O}}_{\mathrm{3}}}_i}.\end{array}$$

This is converted from $g_{{{\mathrm{O}}_{\mathrm{3}}}_i}$ in Eq. (33) by dividing the conductance value in mmol m⁻¹ s⁻¹ by 41 000 (assuming standard temperature (20 °C) and air pressure (1.013×10⁵ Pa)) to give conductance in m s⁻¹.

2.5 Crop biomass, LAI, height, and yield variables

The following section describes how to estimate crop biomass, important canopy characteristics (LAI and crop height (h)), and yield variables from accumulated calculations of Anet_c over the course of the growing season following (Osborne et al., 2015).

2.5.1 Crop biomass (NPP and GPP)

The simulation of crop growth requires an estimate of the net primary productivity (NPP) which is calculated at the end of each day and summed over the growing season. Carbon is assumed to be allocated to five key crop components: root, leaf, stem, harvest, and reserve pools (Osborne et al., 2015). This carbon allocation is ultimately used to simulate leaf area index (LAI), canopy height (h), biomass, harvest index, and yield at the end of each day throughout the growing season.

NPP (kg C m⁻² d⁻¹) is accumulated throughout the day using the JULES-crop approach to model crop growth (Osborne et al., 2015) described in Eq. (35),

$$\begin{array}{}\text{(35)}& \mathrm{NPP}=\mathrm{GPP}-R_{\mathrm{p}},\end{array}$$

where GPP is the gross primary productivity (kg C m⁻² d⁻¹) and R_p is plant respiration divided into maintenance (R_pm) and growth (R_pg) respiration (kg C m⁻² d⁻¹) (Clark et al., 2011), where $R_{\mathrm{p}}=R_{\mathrm{pm}}+R_{\mathrm{pg}}$ and where R_pg is assumed to be a fixed fraction of the NPP as shown in Eq. (36),

$$\begin{array}{}\text{(36)}& R_{\mathrm{pg}}=R_{\mathrm{gcoeff}}(\mathrm{GPP}-R_{\mathrm{pm}}),\end{array}$$

where R_gcoeff is the growth respiration coefficient which was initially set to 0.25 based on the value for all PFTs (i.e. forests and grasses including crops) in Clark et al. (2011). GPP is calculated by Eq. (37):

$$\begin{array}{}\text{(37)}& \mathrm{GPP}={\mathrm{Anet}}_{\mathrm{c}}+f_{\mathrm{PAW}}{R}_{\mathrm{dc}},\end{array}$$

where Anet_c is net canopy photosynthesis (see Eq. 28) and f_PAWR_dc is the soil-moisture modified canopy dark respiration (kg C m⁻² d⁻¹), where $R_{\mathrm{dc}}=V_{\mathrm{cmax},i}\times R_{\mathrm{dcoeff}}$, with R_dcoeff initially assumed to be 0.015 based on Clark et al. (2011); V_cmax,i is the maximum carboxylation efficiency for each canopy layer i, which decreases from the top to bottom of the canopy (see Eq. 30); and f_PAW is calculated in Eq. (7).

Leaf maintenance respiration (R_pm) is assumed to be equivalent to the soil moisture modified canopy dark respiration, while root and stem respiration are assumed to be independent of soil moisture but to have the same dependencies on C content. We assume a fixed relationship between C and N contents of these organs so that R_pm can be estimated by Eq. (38):

$$\begin{array}{}\text{(38)}& R_{\mathrm{pm}}=R_{\mathrm{dc}}\times \left(f_{\mathrm{sw}}+\left({\displaystyle \frac{{\mathrm{C}}_{\mathrm{root}}+{\mathrm{C}}_{\mathrm{stem}}}{{\mathrm{C}}_{\mathrm{leaf}}}}\right)\right).\end{array}$$

The C accumulating as NPP each day is divided into five carbon pools, i.e. root (C_root), leaf (C_leaf), stem (C_stem), reserve (C_resv), and harvest (C_harv) (kg C m⁻² d⁻¹), according to partition coefficients (see Eq. 39) allowing for accumulation of C in these pools over the course of the crop growth period:

$$\begin{array}{}\text{(39)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{dC}\mathit{\_}\mathrm{root}}{\mathrm{d}t}}=p_{\mathrm{root}}\mathrm{NPP},\\ & {\displaystyle \frac{\mathrm{dC}\mathit{\_}\mathrm{leaf}}{\mathrm{d}t}}=p_{\mathrm{leaf}}\mathrm{NPP},\\ & {\displaystyle \frac{\mathrm{dC}\mathit{\_}\mathrm{stem}}{\mathrm{d}t}}=p_{\mathrm{stem}}\mathrm{NPP}\left(\mathrm{1}-\mathit{\tau }\right),\\ & {\displaystyle \frac{\mathrm{dC}\mathit{\_}\mathrm{harv}}{\mathrm{d}t}}=p_{\mathrm{harv}}\mathrm{NPP},\\ & {\displaystyle \frac{{\mathrm{dC}}_{\mathrm{resv}}}{\mathrm{d}t}}=p_{\mathrm{stem}}\mathrm{NPP},\mathit{\tau },\end{array}\end{array}$$

where τ is the fraction of stem C that is partitioned into the reserve pool. p_rootp_leaf, p_stem, and p_harv=1. The partition coefficients are related to the crop development stage (DVI) and hence effective thermal time (TT_eff) since emergence. The partition coefficients are based on Osborne et al. (2015) and provided as a function of DVI using six parameters to continuously describe varying partition coefficients over the duration of the crop growing season. We use the same multinomial logistic as that described in Osborne et al. (2015) to define this function according to Eq. (40):

$$\begin{array}{}\text{(40)}& \begin{array}{rl}& p_{\mathrm{root}}={\displaystyle \frac{e^{{\mathit{\alpha }}_{\mathrm{root}}+\left({\mathit{\beta }}_{\mathrm{root}}\mathrm{DVI}\right)}}{\begin{array}{c}{e}^{{\mathit{\alpha }}_{\mathrm{root}}+\left({\mathit{\beta }}_{\mathrm{root}}\mathrm{DVI}\right)}+e^{{\mathit{\alpha }}_{\mathrm{stem}}+\left({\mathit{\beta }}_{\mathrm{stem}}\mathrm{DVI}\right)}+\\ e^{{\mathit{\alpha }}_{\mathrm{leaf}}+\left({\mathit{\beta }}_{\mathrm{leaf}}\mathrm{DVI}\right)}+\mathrm{1}\end{array}}},\\ & p_{\mathrm{stem}}={\displaystyle \frac{e^{{\mathit{\alpha }}_{\mathrm{stem}}+\left({\mathit{\beta }}_{\mathrm{stem}}\mathrm{DVI}\right)}}{\begin{array}{c}{e}^{{\mathit{\alpha }}_{\mathrm{root}}+\left({\mathit{\beta }}_{\mathrm{root}}\mathrm{DVI}\right)}+e^{{\mathit{\alpha }}_{\mathrm{stem}}+\left({\mathit{\beta }}_{\mathrm{stem}}\mathrm{DVI}\right)}+\\ e^{{\mathit{\alpha }}_{\mathrm{leaf}}+\left({\mathit{\beta }}_{\mathrm{leaf}}\mathrm{DVI}\right)}+\mathrm{1}\end{array}}},\\ & p_{\mathrm{leaf}}={\displaystyle \frac{e^{{\mathit{\alpha }}_{\mathrm{leaf}}+\left({\mathit{\beta }}_{\mathrm{leaf}}\mathrm{DVI}\right)}}{\begin{array}{c}{e}^{{\mathit{\alpha }}_{\mathrm{root}}+\left({\mathit{\beta }}_{\mathrm{root}}\mathrm{DVI}\right)}+e^{{\mathit{\alpha }}_{\mathrm{stem}}+\left({\mathit{\beta }}_{\mathrm{stem}}\mathrm{DVI}\right)}+\\ e^{{\mathit{\alpha }}_{\mathrm{leaf}}+\left({\mathit{\beta }}_{\mathrm{leaf}}\mathrm{DVI}\right)}+\mathrm{1}\end{array}}},\\ & p_{\mathrm{harv}}={\displaystyle \frac{\mathrm{1}}{\begin{array}{c}{e}^{{\mathit{\alpha }}_{\mathrm{root}}+\left({\mathit{\beta }}_{\mathrm{root}}\mathrm{DVI}\right)}+e^{{\mathit{\alpha }}_{\mathrm{stem}}+\left({\mathit{\beta }}_{\mathrm{stem}}\mathrm{DVI}\right)}\\ +e^{{\mathit{\alpha }}_{\mathrm{leaf}}+\left({\mathit{\beta }}_{\mathrm{leaf}}\mathrm{DVI}\right)}+\mathrm{1}\end{array}}},\end{array}\end{array}$$

where DVI is the development index, and α and β are partition parameters. These parameters describe the shape of the thermal-time varying partition coefficient for leaves, roots, and stems.

Once C is no longer partitioned to stems, C from the stem reserve pool will mobilise to the harvest pool at a rate of 10 % per day following Osborne et al. (2015) and described by Eq. (41):

$$\begin{array}{}\text{(41)}& \begin{array}{rl}{C}_{\mathrm{harv}}& =C_{\mathrm{harv}}+\left(\mathrm{0.1}{C}_{\mathrm{resv}}\right)C_{\mathrm{resv}}=\mathrm{0.9}{C}_{\mathrm{resv}}\mathit{\}}\\ & \text{for}{p}_{\mathrm{stem}}<\mathrm{0.01}.\end{array}\end{array}$$

Total leaf C is divided between green leaf C (C_{leaf, green}) and brown leaf carbon (C_{leaf, brown}). Carbon from the C_{leaf, green}  will mobilise to the harvest pool at a rate of 5% per day after Osborne et al. (2015) and to the C_{leaf, brown} at a rate of 24 % per day once f_LS >1 as described in Eq. (42):

$$\begin{array}{}\text{(42)}& \begin{array}{rl}\mathit{\{}{\mathrm{C}}_{\mathrm{harv}}& ={\mathrm{C}}_{\mathrm{harv}}+\left(\mathrm{0.05}{\mathrm{C}}_{\mathrm{leaf},\mathrm{green}}\right){\mathrm{C}}_{\mathrm{leaf},\mathrm{green}}=\mathrm{0.86}\\ & {\mathrm{C}}_{\mathrm{leaf}}{\mathrm{C}}_{\mathrm{leaf}}=\mathrm{0.86}{\mathrm{C}}_{\mathrm{leaf},\mathrm{green}}+\mathrm{0.24}{\mathrm{C}}_{\mathrm{leaf},\mathrm{brown}}\mathit{\}}\\ & \text{for}{f}_{\mathrm{LS}}>\mathrm{1}.\end{array}\end{array}$$

2.5.2 Leaf area index (LAI) and stem height (h)

At the end of each day, the C content of the stem and leaf is used to estimate LAI by Eqs. (43) and (44):

$$\begin{array}{}\text{(43)}& {\displaystyle }\mathrm{LAI}=({\mathrm{C}}_{\mathrm{leaf}}/f_{\mathrm{c}})\times \mathrm{SLA},\text{(44)}& {\displaystyle }\text{where}\mathrm{SLA}=\mathrm{{\rm Y}}(\mathrm{DVI}+\mathrm{0.06}{)}^{\mathit{\delta }}.\end{array}$$

The values Υ and δ were determined by fitting the values to the paired values of DVI and specific leaf area (SLA). The value of f_c is 0.5 (unitless) and denotes the carbon fraction of dry matter.

The amount of C in the stem is used to calculate the crop height h in metres by Eq. (45):

$$\begin{array}{}\text{(45)}& h=k(C_{\mathrm{stem}}/f_{\mathrm{c}}{)}^{\mathit{\lambda }},\end{array}$$

where k and λ were determined by fitting the value C_stem and h.

2.5.3 Yield variables

According to Osborne et al. (2015) yield can be calculated from the C allocated to the harvest pool (C_harv) at the end of the growing season as described in Eq. (46):

$$\begin{array}{}\text{(46)}& {\mathrm{Yield}}_{\mathrm{grain}}={\displaystyle \frac{({\mathrm{C}}_{\mathrm{harv}}\times (\mathrm{1}/f_{\mathrm{c}})\times D_w\times E_{\mathrm{g}})}{\mathrm{1000}}},\end{array}$$

where harvested C is converted to total biomass (using the conversion factor f_c=0.5), i.e. by multiplying the harvested C by $\mathrm{1}/f_{\mathrm{c}}$ and then by 1/0.84 (D_w) to account for the grain moisture content (Mulvaney and Devkota, 2020). C_harv includes both chaff and grain; however, O₃ fumigation experimentalists tend to only include grain when calculating total crop yield at the end of the growing season, so we assume 15 % of the yield is chaff and include a grain-to-ear ratio, E_g, of 0.85. Dividing by 1000 converts yield from kg C m⁻² to g C m⁻², the unit most often used to describe experimental yield results.

Evaluation of the DO₃SE-Crop model uses a variety of growth dry matter (DM) metrics. Some of the most important metrics and their calculations are “straw DM”, which is calculated as the sum of carbon allocated to C_stem, C_leaf, and C_resv; “ear DM”, which is calculated from C_harv excluding the moisture content (D_w) conversion; “grain DM”, which is calculated from C_harv excluding both the moisture content (D_w) conversion and removing the chaff fraction conversion E_g; “above-ground DM”, which is the straw DM plus the ear DM; “below-ground DM”, which is converted from C_root; and “harvest index”, which is the grain DM divided by the above-ground DM. In all cases the f_c conversion factor is used to convert from, for example, g C m⁻² to g DM m⁻².

3 DO₃SE-Crop model calibration

3.1 Xiaoji China experimental dataset

The DO₃SE-Crop model was used to analyse the O₃-FACE (Free Air Concentration Enrichment) experimental data collected in Xiaoji, Jiangdu, Jiangsu Province, China. The wheat crop was grown in fully open-air field conditions for three consecutive growing seasons from 2007 to 2009. The dataset includes four modern cultivars of winter wheat (Triticum aestivum L.) grown under ambient (AA) and elevated (E) O₃, with the elevated treatment being, on average, 25 % above the ambient O₃ concentrations from early March/April to the end of May each year. The four cultivars were Yannong 19 (strong-gluten wheat, hereafter Y19), Yangmai 16 (medium-gluten wheat, hereafter Y16), Yangmai 15 (weak-gluten wheat, hereafter Y15), and Yangfumai 2 (weak-gluten wheat, hereafter Y2) (Zhu et al., 2011).

Soil water availability was sufficient for optimum wheat crop growth, so we assumed there was no soil moisture stress (Feng et al., 2012). Any data gaps were filled following the AgMIP-O₃ gap-filling protocol (see Sect. S4). For large O₃ data gaps (i.e. greater than 2 weeks) occurring outside the O₃ fumigation period, we used scaled WRF-Chem (version 4.2) data for Xiaoji (Conibear et al., 2018) to ensure consistency in model calibration and potential applications across China. The dataset provides grain yield components, including the number of ears per square metre, the number of grains per ear, and the grain dry matter (grain DM, in g m⁻²) (Feng et al., 2011, 2016). Additional physiological datasets (i.e. A_net, V_cmax, J_max, and $g_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ (converted to $g_{{\mathrm{O}}_{\mathrm{3}}}$ as described in Sect. S1b)) are also provided, but only for the year 2008 for all cultivars (Y2, Y19, Y15, and Y16) and for the flag leaf. The 2008 data also include measurements of chlorophyll (in mg m⁻²) which can be used to assess the level of senescence experienced by the leaf (Mariën et al., 2019). Since the year 2008 also showed significant differences in grain DM between AA and E O₃ treatments (a mean relative yield difference of 6.73 for all cultivars; see Table S2b), this year was used to train the DO₃SE-Crop model with other years (i.e. 2007 and 2009) used to test the model.

Further experimental details are provided in Feng et al. (2011, 2016). Table 1 describes the average, minimum, and maximum values for all measured variables required to run the DO₃SE-Crop model collected at the Xiaoji site for each year. Additionally, the M7 (mean 7 h O₃ concentration over the exposure period in ppb) is included for both AA and E O₃ treatments. Measurements were taken at a height of 2 m above the ground surface.

Table 1 Summary of hourly meteorological and ozone concentration ([O₃]) data at Xiaoji.

3.2 DO₃SE-Crop calibration and evaluation

Development and calibration of the DO₃SE-Crop model with the Xiaoji experimental dataset followed three main steps: (i) sensitivity analysis to identify key model parameters to calibrate, (ii) calibration of these key parameters for a single year and both tolerant and sensitive cultivars, and (iii) evaluation of key DO₃SE-Crop model outputs for different years and cultivars from those used in model calibration.

To perform the sensitivity analysis we used the SALib Python library (Iwanga et al., 2022; Herman and Usher, 2017). The analysis requires ranges to be specified for the parameters (identified by an initial manual calibration) that are included in the sensitivity analysis. For physiological parameters, ranges were determined by considering the range of these parameters in the literature. For carbon allocation parameters, the range was identified by considering the maximum and minimum values of these parameters that would result in appropriate dry matter partitioning within the plant. Once the ranges were identified, the sensitivity analysis was run using the extended Fourier amplitude sensitivity analysis, which has been commonly used by other crop modellers to improve their calibrations (Silvestro et al., 2017; Vazquez-Cruz et al., 2014). From the sensitivity analysis outputs (see Fig. S6), the parameters whose variation contributes the most to variations in selected modelling outputs (in this case photosynthetic rate and yield) were identified as the key model outputs for calibration. Using this method we identified the following DO₃SE-Crop parameters as those most important to calibrate: (i) leaf photosynthesis parameters (V_cmax25, J_max25, kN, m, and VPD₀) and (ii) C allocation parameters (α_root, α_leaf, α_stem, Υ , τ) and related dark respiration coefficients (R_dcoeff and R_gcoeff) which were later included in the calibration after identifying issues with overestimated respiration, likely due to the use of parameter values designed for broad plant functional types, which may not be suitable for wheat. O₃ damage module parameters related to senescence (γ₃, γ₄, γ₅, and CLsO3) were not included in the sensitivity analysis, as γ₃ and CLsO3 are already recognised as important for calibration, and γ₄ and γ₅ were introduced in this study to represent the start (SOS) and end (EOS) of senescence, making both essential for calibration. Phenology parameters were also excluded as earlier studies have shown these are relatively straightforward to calibrate using automated methods for a range of environmental conditions (Nguyen et al., 2024). We note that assessing the probability distribution of these ranges would also be useful but consider this outside the scope of the current paper due largely to data limitations.

The DO₃SE-Crop model was then calibrated using the 2008 dataset for the Y2 and Y16 cultivars. The year 2008 was selected since this showed a substantial difference in yield of 208 and 148 g m⁻² between the AA and EO₃ treatments for the Y2 and Y16 cultivars, respectively. These cultivars were chosen since they were identified as the most sensitive (Y2) and tolerant (Y16) cultivars according to the experimental analysis conducted by Feng et al. (2016). See Fig. 5, which shows a diagram representing the calibration process. Calibration of the phenology module used only the Y2 cultivar, AA O₃ treatment data describing the timing of emergence, anthesis, and maturity to calibrate key phenology parameters (T_b, T₀, T_m, VT_min, VT_max, PIV and PID, TT_emr, TT_veg, TT_rep, and T_l). The phenology calibration was automated by computationally applying a genetic algorithm (Wang, 1997), an optimisation technique with gradient decent to find the best parameters. This uses a combination of crossover strategy (selecting parameters randomly from parameter pairings) and mutation strategy (which takes a parameter range and uses incremental step changes) to identify the parameters which give the highest R² and lowest root mean square error (RMSE) when compared with observations of the timing (day of year) of anthesis and maturity.

Calibration of the leaf physiology, canopy C allocation, and O₃ damage DO₃SE-Crop modules was performed manually. This required that an initial value and range be defined for each parameter, which were defined from a combination of observations from the Xiaoji experimental dataset as well as values taken from the literature (see Tables A1 and A2 of the Appendix A for details). The model was manually calibrated until certain conditions were satisfied, as explained below. Calibration of the leaf physiology parameters (V_cmax, J_max, kN, m, and VPD₀) was performed only for the Y2 cultivar, AA O₃ treatment whilst keeping all other parameters fixed. This calibration aimed to achieve a maximum A_net value of 30 µmol CO₂ m⁻² s⁻¹ and a $g_{{\mathrm{O}}_{\mathrm{3}}}$ value of 350 mmol O₃ m⁻² PLA s⁻¹, consistent with the maximum values observed in the Xiaoji dataset (Zhu et al., 2011). We calibrated V_cmax and J_max as measurements are only provided for Y2 and Y16 cultivars and only for certain points during the growth period, and we know that V_cmax and J_max can vary seasonally.

Calibration of the C allocation parameters (α_root, α_leaf, α_stem, Υ, τ) and related dark respiration coefficients (R_dcoeff and R_gcoeff) was also performed keeping all other parameters fixed. This calibration aimed to achieve the following criteria: a stem dry matter to leaf dry matter ratio (R_SL) of approximately 2:1 (Huang et al., 2022), relative growth of different plant parts (i.e. leaves, stem, roots, grain) consistent with profiles found in the literature (Osborne et al., 2015; Penning de Vries, 1989), a modelled grain DM within ±30 % of the observed, an above-ground DM value of between 1200–1600 g m⁻², an LAI value between 4–7 m² m⁻², and an R_d value of between 30 % and 60 % of A_net (Amthor et al., 2019). We calibrated C allocation parameters as the JULES-crop model calibration has only been performed for broad, global-scale application for wheat (Osborne et al., 2015) and therefore requires further calibration for application under Chinese conditions. Further, the observed dataset does not provide any information with regards to the change in carbon allocation parameters due to ozone. The C allocation parameters were only calibrated for ambient ozone conditions, and we only investigate the effect of ozone on C assimilation (not C allocation).

Finally, calibration of the O₃ parameters (γ₃, γ₄, and γ₅) was performed using 2008 data for the Y2 and Y16 cultivars whilst again keeping the other parameters fixed. Calibration was targeted so that the difference in grain DM between ambient and elevated O₃ treatments was as close as possible to ±10 % of the observed.

The manual calibration process consisted of three stages as explained above, as well as comparisons with established information on wheat growth from the literature. By reducing the number of parameters involved in the calibration, the chance of equifinality (multiple combinations of parameters yielding similar results) was minimised (Beven, 2006). The parameters identified by the sensitivity analysis were varied within realistic ranges to obtain a parameterisation that closely approximates wheat physiological processes. Multiple parameterisations were tested to avoid convergence on local minima in R² and RMSE. While further fine-tuning of the parameter ranges could potentially improve yield prediction, it might also disrupt simulations of other key plant processes, such as carbon allocation or photosynthesis. The calibration approach balances the need for accurate output simulation with the physiological realism required for wheat growth under the conditions of this study. Though it is difficult to claim that the absolute optimal parameter set has been achieved, this limitation is common to any model calibration (Wallach, 2011). The current parameterisation represents a physiologically realistic simulation of wheat growth under the conditions of the present study using a robust calibration method.

Evaluation of the DO₃SE-Crop model was conducted using Xiaoji data for 2007 and 2009 for all cultivars and 2008 data for Y19 and Y16 cultivars. This evaluation tested the ability of the calibrated DO₃SE-Crop model to simulate grain DM using R² and RMSE statistical tests.

4 Results

We first examine the model's ability to simulate the key phenological development stages since this is key to simulating the variation in C allocation to different plant parts over the course of the growing season and hence how O₃ exposure will influence growth and yield which is determined by the timing and length of the grain-filling period. We also explore how DO₃SE-Crop simulates within-canopy [O₃] profiles to understand which layers of the canopy are most important in determining O₃ response. We then examine the ability of the model to simulate leaf-level physiology and C allocation to the different parts of the crop. Lastly, the impact of both instantaneous and long-term O₃ damage on the crop's final grain DM is evaluated for different cultivars and years.

i. Crop phenology

The Xiaoji dataset provides sowing and harvest dates for all cultivars for each year but only provides the date of the timing of anthesis for the years 2008 and 2009 for all cultivars. We assume that DVI =1 is equivalent to the start of anthesis and that this occurs 4–5 d after flag leaf emergence as shown in Fig. 2. We determine the influence of O₃ on the start and end of senescence (SOS and EOS) using the breakpoint method (described in Pande et al., 2024) to assess significant changes in the chlorophyll values that indicate senescence onset and rate of change for the quantification of tl_ep and tl_se. This method is applied to chlorophyll data collected in 2008 under both AA and E O₃ treatments for the Y2 cultivar. We then assume that these key phenology parameters (i.e. TT_emr, TT_veg, TT_rep, tl_ep, and tl_se) are consistent across cultivars and years. Our results in Fig. 3 suggests this is a reasonable assumption; however, we appreciate that assuming these phenology parameters will work for a wider variety of cultivar types (e.g. early or late sown and/or maturing) and years with rather different meteorological conditions needs to be done with caution.

Figure S1 shows the modelled vs. observed timing of anthesis and harvest for the training dataset. Figure 3 shows the same for the test dataset. For the test dataset there is a variation of 2 to 4 and 1 to 6 d for the modelled anthesis and maturity in relation to observed anthesis and maturity, respectively, with observed phenology tending to be a little later than modelled. The T_l ranges between 1325 and 1478 °C days for the 3 years, with crop sowing occurring between 315 and 324 d of year and harvests occurring between 135 and 151 d of year (of the following year). The number of days from the modelled crop sowing to harvest was between 181 and 191 for the 3 years, compared to 198 and 201 for the observations.

Figure 3 Modelled vs. observed phenological stages provided as day of year (DOY) for the test dataset (i.e. excluding the year 2008 for the Y2 cultivar).

Figure 4 Comparison of daily maxima seasonal profiles of DO₃SE-Crop modelled canopy leaf vs. observed flag leaf data for (a) AA O₃ treatment A_net, (b) AA O₃ treatment $g_{{\mathrm{O}}_{\mathrm{3}}}$, (c) E O₃ treatment A_net, and (d) E O₃ treatment $g_{{\mathrm{O}}_{\mathrm{3}}}$ for the period from the anthesis (i.e.  TT_rep) for the year 2008 and the Y16 cultivar. The left (solid blue line) and right (solid red line) represent the segment fits to the normalised chlorophyll content values for application of the breakpoint method to define the SOS (start of senescence) shown as the dashed black line. The green scatter solid dots, along with their standard measurement error, represent the normalised observed chlorophyll content values (see Fig. 7 for further details).

ii. Leaf physiology variables (A_net, $g_{{\mathrm{O}}_{\mathrm{3}}}$)

The DO₃SE-Crop model was able to simulate the seasonal A_net and $g_{{\mathrm{O}}_{\mathrm{3}}}$ with values ranging from 0 to 27 µmol CO₂ m⁻² s⁻¹ and 10 to 351 mmol O₃ m⁻² s⁻¹ for A_net and $g_{{\mathrm{O}}_{\mathrm{3}}}$, respectively, over the course of the growing season (see Fig. 4). The simulated daily maximum values of modelled $g_{{\mathrm{O}}_{\mathrm{3}}}$, at 351 mmol O₃ m⁻² s⁻¹, were within the range of the observed value of 340 mmol O₃ m⁻² s⁻¹. Similarly, the modelled daily maximum A_net is 27 µmol CO₂ m⁻² s⁻¹ compared to observed value of 28 µmol CO₂ m⁻² s⁻¹ for the period between anthesis and 10 d before maturity for the year 2008, for the Y16 cultivar (similar results were obtained for the Y2 cultivar; see Fig. S5). In Fig. 4a and b, the steep decline in modelled A_net and $g_{{\mathrm{O}}_{\mathrm{3}}}$ is not seen in the observed dataset. This discrepancy may occur since the simulated A_net and $g_{{\mathrm{O}}_{\mathrm{3}}}$ values represent sunlit parts of the upper canopy which comprise both green and senesced leaf material. In contrast, observed A_net and $g_{{\mathrm{O}}_{\mathrm{3}}}$ values are measured specifically on the flag leaf and most likely only for the green parts of the leaf, since the LI-6400 photosynthesis system mounted with a 6400–40 leaf chamber fluorometer (used to measure A_net and $g_{{\mathrm{O}}_{\mathrm{3}}}$ in the Xiaoji experiment, Feng et al., 2016) will not provide values for senesced leaf material. See also Fig. 4, which combines A_net and $g_{{\mathrm{O}}_{\mathrm{3}}}$ with observed normalised chlorophyll content and clearly shows the leaf is senescing as predicted by the model. However, the observed decline in chlorophyll values closely matches the decrease in modelled A_net and $g_{{\mathrm{O}}_{\mathrm{3}}}$, with the model accurately capturing the timing of the earlier onset of senescence, which occurred 0–3 d earlier in the AA and EO₃ treatments. It is useful to note that the calibrated V_cmax and J_max values match the observed values within ±2 µmol CO₂ m² s⁻¹.

Figure 5 Plot showing variation in key O₃ deposition terms as daily maxima by canopy layer (NB i=4 is the top canopy layer, n=4) (a) O₃ concentration at the top of each layer, (b) leaf boundary layer resistance by canopy layer ($rb{,}_{{\mathrm{O}}_{\mathrm{3}}}$), (c) PAR for the sunlit LAI component of each layer (PAR_sun), and (d) leaf-level stomatal conductance to O₃ ($g_{{\mathrm{O}}_{\mathrm{3}}}$) for the period from anthesis (i.e. TT_rep) for the Y16 cultivar and for the E O₃ treatment in 2008.

iii. Within-canopy variation in O₃ and physiology

An important determinant of O₃ deposition and damage is stomatal O₃ deposition (our $g{{\mathrm{O}}_{\mathrm{3}}}_{\mathrm{c}}$), which is a function of within-canopy transfer of O₃ and stomatal and non-stomatal deposition. The multi-layer aspect of the DO₃SE-Crop model allows within-canopy stomatal and non-stomatal O₃ deposition to be simulated. Figure 5 shows the variation in key variables that determine total and stomatal O₃ canopy deposition across four canopy layers as a midday average over the course of the tl_ep period of the flag leaf, for the year 2008 and the Y16 cultivar.

Table 2 Simulations of percent of grain DM loss that compare ambient (AA) and elevated (E O₃) treatments with a pre-industrial O₃ scenario divided between grain DM losses caused by the instantaneous effect O₃ on photosynthesis and the long-term O₃ effect on senescence. The effects of both damage O₃ mechanisms acting together are also shown.

Figure 5a shows a decrease of within-canopy O₃ concentration from highs of around 140 ppb to values within the range of 10 to 50 ppb between the top of the canopy and bottom canopy layer; the penetration of O₃ into the canopy increases over time as the canopy senescence and O₃ uptake is reduced. Similarly, PAR_sun is reduced from maximum values of around 200 W m⁻² at the top of the canopy to values of around 100 W m⁻² in the lower canopy layers even on sunny days (see Fig. 5c). The leaf $r{b}_{{\mathrm{O}}_{\mathrm{3}}}$ (Fig. 5b) increases with canopy depth with resistances in the region of approximately 50 s m⁻¹ at the top of the canopy to values of around 600 s m⁻¹ at the bottom of the canopy; this will limit stomatal O₃ uptake in the lower canopy layers; finally, these factors combine to influence canopy level $g_{{\mathrm{O}}_{\mathrm{3}}}$ (Fig. 5d), which is reduced from values of around 350 at the top of the canopy to 20 nmol O₃ m⁻² s⁻¹ at the bottom of the canopy layer; these differences in leaf $r{b}_{{\mathrm{O}}_{\mathrm{3}}}$ and $g_{{\mathrm{O}}_{\mathrm{3}}}$ are reduced with the onset of senescence. This analysis shows the importance of interplay between these different factors for an accurate whole-canopy estimate of O₃ deposition.

iv. Crop development, biomass, and yield

The dry matter dynamics of the different parts of the crop are shown in Fig. 6. The modelled grain DM value of 851 g m⁻² was reasonably close to the observed value of 888 g m⁻². The stem to leaf dry matter ratio (R_SL) is 2.1:1 and therefore in the range provided in the literature (Huang et al., 2022). The above-ground biomass values of 1510 g m⁻² also match reasonably well against the 1200 to 1600 g m⁻² range described in the literature (Huang et al., 2022; Liu et al., 2022). Further, the partition fraction profiles are consistent with those of Osborne et al. (2015) as shown in Fig. 6a, with the main differences being that the modelled stem and root partition profiles are somewhat higher and lower, respectively. The JULES model comparison is provided for illustrative purposes only (i.e. this model has not been calibrated with the Xiaoji data but rather is a parameterisation suggested for global application).

Figure 6 Seasonal profiles (i.e. plotted against DVI) of carbon allocation variables for the Xiaoji calibrated DO₃SE-Crop model (i.e. AA O₃ treatment, year 2008, and Y16 cultivar) with panel (a) showing the partition fractions of the daily accumulated NPP partitioned to roots, stems, leaves, and grains for the Xiaoji calibrated DO₃SE-Crop model (solid lines) vs. the JULES-crop model (dashed line) calibrated for global application after Osborne et al. (2015) and panel (b) showing the DM (in g m⁻²) of daily accumulated NPP partitioned to roots, stems, leaves, and grains, with the observed final grain DM for Y16 cultivar in 2008 also shown (solid black dot with 5 % error).

v. O₃-induced yield loss difference between tolerant and sensitive cultivars: instantaneous and long-term senescence impact

The grain DM is assumed to be damaged by both the instantaneous impact of O₃ (Farage et al., 1991) on photosynthesis and a longer-term O₃ effect that can lead to enhanced senescence (Y. Feng et al., 2022). To explore which of these damage mechanisms is most important, we calculated the difference in the grain DM caused by carbon assimilation for the AA and E O₃ treatments as compared to a simulated very low O₃ treatment representing pre-industrial conditions (for which C_z O₃ concentration did not exceed 15 ppb) for the tolerant (Y16) and sensitive (Y2) cultivar for each of the 3 years (see Table 2). We found a negligible effect of O₃ (0 % to 0.2 %) on grain DM due to the instantaneous effect of O₃ on photosynthesis, which could perhaps be partly due to the to the crops ability to recover photosynthetic capacity overnight, compared to a highly significant (9.85 % to 31.13 %) impact due to the long-term O₃ effect on carbon assimilation via the enhancement of senescence on final grain DM. Table S3 shows the observed percent of grain DM loss compared to a modelled pre-industrial O₃ scenario due to the combination of instantaneous and long-term ozone effect.

vi. Senescence

The breakpoint method (Mariën et al., 2019) was used to determine the onset (SOS) and end (EOS) of senescence and maturity, respectively, using the chlorophyll data which were available for the year 2008 and the Y16 and Y2 cultivars. Results in Figs. 7 and S4 show that the E O₃ treatment for cultivars Y16 and Y2 brought forward the SOS by 3 and 5 d (see Fig. 7), respectively, and EOS by 6 and 9 d (see Fig. S4), respectively. Figure 7 also shows the f_LS profile which denotes the DO₃SE-Crop models' accumulated stomatal O₃ flux effect on senescence; it is clear that f_LS is able to simulate the change in normalised chlorophyll content reasonably well. The slope of the ambient f_LS is already steep since the ambient treatment already has rather high O₃ levels as is now made clear in Table 1 with a value of 47 ppb. According to the M7 wheat dose–response relationship this would result in a yield loss of ∼5 %.

Figure 7 Profiles of O₃-induced leaf senescence for the Y16 cultivar for the (a) AA O₃ treatment and (b) E O₃ treatment. The timing of the SOS (solid black line) and EOS (dashed black line) was determined by applying the breakpoint method to the chlorophyll data and is shown in relation to the f_LS simulations of senescence (solid yellow line). The observed normalised chlorophyll content data, shown as filled green symbols, include error bars representing the standard deviation of the measurements.

Figure 8 Box plots (crosses: 0.01 and 0.99 percentiles; box: 0.25 quartile, median, and 0.75 quartile; triangle: mean) of simulated and observed wheat grain DM for the tolerant (Y15 and Y16) and sensitive (Y2 and Y19) cultivars under (a) AA and (b) E O₃ treatment for the years 2007, 2008, and 2009; these data include the whole dataset.

Figure 9 A scatter plot showing modelled vs. observed grain DM (in g m⁻²) for the AA and E O₃ treatments for all four cultivars and 3 years of the Xiaoji dataset; these data include those used for evaluation.

vii. Grain DM simulations across years and between cultivars

Figure 8 shows a box plot of the modelled vs. observed grain DM for the sensitive (Y2, Y19) and tolerant (Y15, Y16) cultivars for each O₃ treatment (AA and E) for the years 2007, 2008, and 2009 (i.e. all data). Given the variability in the experimental data, the model simulates the difference in grain DM between the AA and E O₃ treatments reasonably well with a simulated reduction in grain DM of 29 to 131 g m⁻² compared with observed values of 81 to 165 g m⁻² for the tolerant cultivars and 49 to 196 g m⁻² compared with observed values of 54 to 293 g m⁻² for the sensitive cultivars, respectively. The most notable difference is that there is a larger range in the simulated grain DM losses of the modelled sensitive cultivars, though the simulated mean value for absolute grain DM suggests a more conservative influence of O₃, with yields at 610 g m⁻² vs. observed average yields of 590 g m⁻².

Finally, Fig. 9 shows the relationship between modelled vs. observed grain DM (in g m⁻²) as a scatter plot; a linear regression through these data gives an R² value of 0.68 and RMSE of 76 g m⁻², showing the model is able to simulate with reasonable accuracy the differences in absolute yield for different cultivars and for different years. There are some instances of both underestimation and overestimation; however, the deviations from the 1:1 line are not excessively large. These models test results compare with an R² of 0.92 (n=4) and an RMSE of 25.49 g m² for the training dataset (Y2 and Y16 cultivar and year 2008; see Fig. S3); the stronger agreement between observed and modelled training dataset, as well as the reasonable agreement for the entire dataset, would suggest the model is not over-fitted. We find that we tend to underestimate the O₃-induced relative yield loss (RYL) by between −2.76 and 15.34 % (observed minus modelled RYL) across all years and cultivars.

5 Discussion

The DO₃SE-Crop model was found capable of simulating O₃ damage to grain yield for O₃-FACE conditions at the experimental site in Xiaoji, China, with a good degree of accuracy. Simulated relative yield losses (RYLs) between AA and E O₃ treatments for all years ranged from 11 % to 14 % and from 13 % to 19 % for tolerant and sensitive cultivars, respectively; these tend to be lower (particularly for the more extreme O₃-induced yield losses of the sensitive cultivars) than the observed values of 13 % to 20 % and 10 % to 35 %. Overall, simulations of tolerant and sensitive cultivars underestimated RYLs by 4 % and 7 %, respectively, on average across years and cultivars (see data in Sect. S6). This would suggest that O₃-induced yield losses can be more reliably modelled for tolerant cultivars, possibly because additional processes causing O₃-induced yield losses in sensitive cultivars are not captured. Such processes might include the effect of O₃ on the allocation of carbon to different plant parts (Feng et al., 2008) or O₃ inducing additional respiratory costs via the upregulation of defence mechanisms (Biswas et al., 2008). The model was also able to simulate absolute grain DM reasonably well. Under AA O₃ levels, grain DM simulated for all years and cultivars was between 616 and 851 g m⁻² compared to observations of between 537 and 982 g m⁻². There is a tendency to overestimate grain DM under ambient conditions and underestimate grain DM under elevated O₃, which is reflected in the RYL values.

Overall, the DO₃SE-Crop model simulation results compare favourably to results made by the MCWLA-Wheat model (Tao et al., 2017), which was also calibrated for the Xiaoji experimental conditions but without distinction between tolerant and sensitive varieties. MCWLA-Wheat simulations of absolute yield varied from ∼5700 to 9000 kg ha⁻¹ (compared to ∼5700 to 9800 kg ha⁻¹) for ambient and from ∼4800 to 8000 kg ha⁻¹ (compared to ∼5200 to 8000 kg ha⁻¹) for elevated O₃ treatments. A mean relative yield loss of 14 % was simulated by the model.

It is useful to set these site-specific estimates of O₃-induced yield losses in the context of yield losses estimated using more traditional, concentration-based O₃ risk assessment methods. A seminal paper by Z. Feng et al. (2022) estimated mean relative yield losses across East Asia due to ambient O₃ concentrations at 33 % (with a mean range of 28 % to 37 %) according to a mean monitored O₃ concentrations of 30.9 ppm h expressed as AOT40 (6-month accumulated daytime O₃ concentration above a threshold of 40 ppb). The mean difference in AOT40 (accumulated over only 75 d) between the AA and E O₃ treatments at Xiaoji across all years was 7.8 ppm h, giving a mean relative yield loss of approximately 10 % to 20 % depending on year and cultivar. As such, our modelled results in terms of RYLs between AA and E O₃ treatments are consistent with these broader results for East Asia.

Crop phenology plays a crucial role in determining the timing of the important O₃ exposure period (i.e. from anthesis to maturity) and hence O₃ damage. Evaluation of the DO₃SE-Crop phenology model shows the model is able to accurately simulate crop phenology for the 3 years at Xiaoji (R²=0.95 and RMSE =2.5; see Fig. 3). Estimating the correct timing of anthesis is crucial since the period from anthesis to crop maturity is the O₃-sensitive period. During this period, accumulated stomatal O₃ flux (acc_fst) will contribute to early and enhanced senescence once the critical threshold (CLsO3) is exceeded. This period also coincides with carbon accumulation in the grain (Kohut et al., 1987; Feng et al., 2008), which may be limited by O₃-induced early onset or enhanced senescence. The DO₃SE-Crop model was developed to accommodate the full range of effects of O₃ on senescence with revised functions, similar to those first developed by Ewert and Porter (2000), and was able to modify both the O₃-induced onset of senescence and the O₃ effect on maturity. This is important since experimental evidence has shown that O₃ can bring forward the maturity date; for example, the flag leaf was found to have senesced 25 d earlier in a high-O₃ treatment, compared to a charcoal-filtered treatment (Grandjean and Fuhrer, 1989; Gelang et al., 2000). O₃ was also found to cause differences in the time to maturity of the flag leaf, with Shi et al. (2009) reporting that maturity was brought forward by 8 d under an elevated O₃ treatment (50 % higher than ambient). Currently, other crop models with O₃ damage functions (e.g. MCWLA-Wheat (Tao et al., 2017) and LINTULCC-2 (Y. Feng et al., 2022)) are only able to bring the O₃-induced onset of senescence earlier.

The DO₃SE-Crop model is also able to simulate differential O₃ uptake in each canopy layer. Figure 5 shows that the majority of stomatal O₃ uptake occurs in the sunlit layers of the upper canopy. Similar results were found in an experimental study on a productive grassland in Switzerland; Jaggi et al. (2006) found that different levels of O₃ exposure to canopy components predominantly located in the upper and lower parts of the canopy support a multi-layer approach to modelling O₃ uptake. Therefore, the focus on the upper canopy by flux-based O₃ metrics (e.g. the phytotoxic ozone dose POD_y; LRTAP, 2017) seems rational in the absence of multi-layer modelling. Crop models such as LINTULCC-2 (Y. Feng et al., 2022) also focus on estimating stomatal O₃ uptake at the top of the canopy to estimate O₃-induced yield losses. For wheat, such an approach is further supported by the fact that the upper canopy layers consist of the flag leaf, which plays a crucial role in photosynthesis and grain filling (Pleijel et al., 2007). The multi-layer functionality of the DO₃SE-Crop model may, however, become more useful when considering crops that partition assimilated carbon to harvest organs earlier in their growing season such as potato (Okrah et al. 2023).

Our results show that the DO₃SE-Crop model was able to estimate the seasonal course of leaf A_net and $g_{{\mathrm{O}}_{\mathrm{3}}}$ daily maxima observed at the Xiaoji site (see Fig. 4a) and when compared to other literature describing leaf physiological variables (Guan et al., 2015; Li et al., 2022). This suggests the coupled A_net–g_sto model is working for Chinese conditions (having previously been applied to and evaluated for European O₃ experimental conditions – see Pande et al., 2024). The leaf physiology parameters used in this study (i.e. for Asian conditions and cultivars) are higher than parameters for European studies. For Europe, V_cmax values of between 60 and 90 µmol CO₂ m⁻² s⁻¹ were found in the literature (Y. Feng et al., 2022; Pande et al., 2024; Van Oijen and Ewert, 1999) compared to the observed mean maximum value of 137 µmol CO₂ m⁻² s⁻¹ at Xiaoji which was used in this study. Similarly, European J_max values ranged from 160 to 180 µmol CO₂ m⁻² s⁻¹ (Feng et al., 2021; Pande et al., 2024; Van Oijen and Ewert, 1999) compared to the observed Xiaoji mean maximum value of 228 µmol CO₂ m⁻² s⁻¹. Even though these leaf physiology parameters are higher, absolute yields for these Chinese cultivars are consistent with those found under European conditions. This most likely reflects the importance of other environmental conditions (e.g. high vapour pressure deficits) limiting leaf carbon assimilation. Moreover, the complex interactions between O₃ exposure and the plants' physiological responses also play a crucial role. Ozone significantly affected antioxidative enzymes, thereby limiting overall photosynthetic efficiency and yield, particularly in O₃-sensitive cultivars, despite their ability to maintain high carboxylation capacity.

Ensuring the seasonal variation in carbon allocation to the different components of the crop (i.e. roots, stem, leaves, and harvest organs) is essential for the simulation of crop growth and yield. There are limited data in the literature that provide these variables, so we compare our results to the carbon allocation profiles described for wheat provided in the original JULES-crop model description, recognising this is intended for wheat grown globally. The DO₃SE-Crop model carbon allocation to the stem and roots is comparatively higher than what was simulated by JULES-crop (Osborne et al., 2015; see Fig. 6a). However, we can justify the carbon allocation coefficients used for Xiaoji since the DO₃SE-Crop model was able to distribute carbon to different plant components to produce a well-proportioned plant over the course of the growing season; this was determined by the calibration to a number of key crop variables (i.e. ratios of plant respiration, LAI, stem to leaf dry matter, above-ground components, and grain dry matter). Importantly, when applied to the test dataset (i.e. excluding 2008 data for the Y2 and Y16 cultivar), the model was found to simulate the grain dry matter under ambient and elevated O₃ treatments to within 7.9 %–8.7 % of the observed values (R²=0.68, 76 g m⁻² see Fig. 9).

The DO₃SE-Crop model, similar to other crop models with O₃ damage functions (i.e. MCWLA-Wheat (Tao et al., 2017), LINTULCC-2 (Y. Feng et al., 2022), and WOFOST (Nguyen et al., 2024)) has the capacity to simulate both the instantaneous and long-term O₃ impact on wheat grain yield. The instantaneous O₃ effect on photosynthesis may cause leaf cell damage and decrease the supply of carbohydrate precursors, which can significantly decrease $g_{{\mathrm{O}}_{\mathrm{3}}}$, V_cmax, and leaf chlorophyll content (Farage et al., 1991). Elevated O₃ also leads to generation of reactive oxygen species (ROS) in plant cells, which can cause oxidative damage to various cellular components. Rubisco, the enzyme responsible for carbon fixation in the photosynthetic process, can be particularly susceptible to this damage, leading to a reduced carboxylation rate (V_cmax). Such an O₃ effect on V_cmax reduces net photosynthesis and can also induce early senescence, shortening the grain-filling period (Triboi and Triboi-Blondel, 2002).

Results from the DO₃SE-Crop model found a larger impact on yield due to the long-term O₃ impact causing relative yield loss of between 10 % and 31 % compared to only between 0 % and 0.2 % resulting from the instantaneous O₃ impact on photosynthesis. Previous studies have also found that the long-term O₃ effect has a larger impact on yield compared to the instantaneous effect of O₃ on photosynthesis (Emberson et al., 2018; Brewster et al., 2024). Senescence is an age-dependent process of degradation and degeneration that allows nutrients to be re-distributed to different plant organs (Lim et al., 2007). Under O₃ stress, this process is often found to occur earlier and more rapidly in leaves as well as at the whole-plant or crop canopy scale (Brewster et al., 2024). The causes of this early and accelerated senescence are not completely understood but may be related to O₃-induced enhanced expression of many genes involved in natural senescence (Miller et al., 1999). Elevated O₃ was also found to inhibit sugar export from leaves (Yadav et al., 2020; Feng et al., 2024), which could trigger early onset of leaf senescence.

The DO₃SE-Crop model accounts for the impact of O₃ on the rubisco enzyme by incorporating modified (Ewert and Porter, 2000) functions for instantaneous and long-term O₃ impact on V_cmax as an important parameter used to characterise the crop photosynthetic capacity (Ewert and Porter, 2000; Osborne et al., 2019). The DO₃SE-Crop model assumes that the O₃ will only accumulate on exceedance of a stomatal O₃ flux threshold of 6 nmol O₃ m⁻² s⁻¹. The long-term O₃ impact mechanism of the DO₃SE-Crop model simulated the effect of senescence on V_cmax reasonably well as evidenced by the reduction in leaf chlorophyll content. We used the breakpoint method (Yang et al., 2016; Mariën et al., 2019) to estimate the SOS and EOS using the day of the year and measured chlorophyll content (Figs. 7 and S4). It is crucial to accurately model the timing of SOS and EOS correctly as this determines the O₃ effect on the duration of the grain-filling period and hence the difference in yield loss due to different O₃ treatments. For example, we modelled a difference of 3 to 5 d in SOS and 6 to 9 d in EOS on average across years for the sensitive and tolerant cultivar, respectively.

China's wheat breeding programme has seen more than 1850 varieties used across China between the 1920s and 2014, leading to increased yields from less than 1 to more than 5 t ha⁻¹ (Qin et al., 2015). Here, albeit with an extremely limited dataset, we parameterise the DO₃SE-Crop model for tolerant and sensitive wheat crop cultivars, since many experimental studies have shown that the response of different cultivars to O₃ stress differs (Biswas et al., 2008). Based on the available data, the model seemed able to capture the difference in grain dry matter between these different cultivar groups across different years reasonably well when compared to the observed dataset (R²=0.68; see Fig. 8). Such a cultivar sensitivity-based parameterisation can provide additional information on the certainty of regional yield loss estimates given the large number of wheat varieties grown across China. However, when applying the model to a broader region, it would be advisable to calibrate phenology for different agro-ecological zones as the temperature changes across China, impacting the duration of the key phenological stages such as anthesis and maturity (Luo et al., 2021). Additionally, carbon allocation parameters may need adjustment, as studies have shown changes in dry matter content across different agro-ecological zones (Hussain and Bangash, 2017).

6 Conclusions

We have shown that the newly developed DO₃SE-Crop model can be calibrated for O₃-tolerant and sensitive wheat varieties for O₃-FACE site conditions at Xiaoji in China. The model can simulate crop phenology, leaf physiology, crop growth, and yield reasonably well across different years. The model is also able to simulate the effect of O₃ stress on grain yield distinguishing the extent of O₃ damage resulting from the same O₃ treatment on cultivars with differing O₃ sensitivities. The DO₃SE-Crop model also has the advantage of simulating O₃ transfer and deposition dynamics within the wheat crop canopy, which could in the future improve our understanding of whole-canopy O₃ effects for crops with different carbon allocation profiles. The ability of the model to estimate relative yield losses across years also suggests the model is fit for purpose to assess the effects of O₃ under a variety of climate variable and O₃ concentration conditions.

Appendix A

Table A1 DO₃SE-Crop variables.

Table A2 DO₃SE-Crop parameters for wheat. DO₃SE-Crop parameters for wheat, including default values. The table includes associated ranges only for parameters that require calibration under varying environmental conditions.