Ozone dose-response relationships for wheat can be derived using photosynthetic-based stomatal conductance models

Ground-level ozone (O 3 ) pollution occurs across many important agricultural regions in Europe, North America, and Asia, negatively impacting O 3 -sensitive crops such as wheat. Risk assessment methods to quantify the magnitude and spatial extent of O 3 pollution have often used dose-response relationships. In Europe, the dose metrics used in these relationships have evolved from concentration-to flux-based metrics since stomatal O 3 flux has been found to correlate better with yield losses. Estimates of stomatal conductance ( g sto ) have to date used an empirical multiplicative model. However, other more mechanistic approaches are available, namely the coupled photosynthetic-stomatal conductance ( A net g sto ) model. This study used a European O 3 OTC and solardome fumigation experimental dataset (comprising 6 cultivars, 4 countries and 14 years) to develop a new flux-based dose-response relationship for wheat yield using the mechanistic A net g sto model (cid:0) A net g sto mech ). The A net g sto mech model marginally improved the regression of the dose-response relationship (R 2 = 0.74) when compared to the flux-response models derived from empirical g sto models. In addition, the A net g sto mech model was somewhat better at predicting the effect of high O 3 concentrations on diurnal and seasonal profiles of g sto and A net . It was also better able to simulate changes of up to 7 and 12 days, respectively, in the start (SOS) and end (EOS) of senescence, an important determinant of yield loss, over a range of O 3 treatments. We conclude that A net g sto mech model can be used to derive robust flux-response relationships.


Introduction
Empirical evidence from Europe, North America and Asia shows that O 3 is causing a range of impacts on staple crops such as wheat (Hansen et al., 2019;Feng et al., 2022;Büker et al., 2015).These impacts include altered stomatal conductance (g sto ) (Danielsson et al., 2003;Ghosh et al., 2020), reduced photosynthesis (A net ) (Ojanperä et al., 1998) and early and enhanced leaf senescence (Osborne et al., 2019;Gelang et al., 2000).Effects on leaf senescence can lead to a reduction in A net and g sto and a shorter grain-filling period (Gelang et al., 2000) thus decreasing yield (Pleijel et al., 2022) and biomass (Feng et al., 2021).Experimental meta-analyses have found that wheat yield losses can range from 3 to 50 % when O 3 concentrations (described as a 7hr daylight mean over the growing season) range from 5 to 115 ppb (Mills et al., 2018).Risk assessments performed on application of dose-response relationships derived from such experimental data (Pleijel et al., 2007) estimate O 3 induced yield losses of between 12 and 15 % globally, causing production losses of approximately 85 million tonnes (Mills et al., 2018).These losses in productivity are a cause for concern, given the importance of wheat as a staple crop for approximately 35 % of the global population (Grote et al., 2021) and that the annual consumption of wheat worldwide is approximately 791 million tonnes (United States Department of Agriculture, 2023).Evidence also suggests that the threat from O 3 pollution will continue into the future.Background O 3 concentrations have remained high over agriculturally important regions (Feng et al., 2019;Arnold et al., 2021;Boleti et al., 2020;Sicard et al., 2021) across Europe (Rega et al., 2020) and both background and peak O 3 concentrations are increasing in the Indo-Gangetic plains in south Asia (Shah et al., 2019), and the North China Plain in East Asia (Liu et al., 2016).To estimate the threat from O 3 pollution, risk assessment modelling methods have been developed to assess the current and future effects of O 3 on crop growth and yield at national, regional, and global scales (Emberson et al., 2018).These methods often use experimental O 3 filtration/fumigation data to derive dose-response relationships and hence require the identification of a suitable dose metric capable of predicting O 3 damage (i.e., yield loss for crops).Metrics would ideally be able to incorporate the effects of species and cultivar as well as management practices (e.g.irrigation) that are known to alter sensitivity to O 3 pollution (Mills et al., 2018;Anav et al., 2016;Osborne et al., 2019).Metrics have evolved over the past decade moving from concentrationto flux-based indices (Grulke and Heath, 2019;Pleijel et al., 2007;Mills et al., 2018) with the flux-based approach allowing O 3 concentrations to be decoupled from O 3 exposure when conditions (e.g., high atmospheric or soil water deficits) limit stomatal O 3 uptake (Emberson et al., 2018;Tai et al., 2021).This capability of the flux-based approach has been shown to give more reliable estimates of the spatial extent of O 3 damage (Mills et al., 2011).
Consequently, the stomatal O 3 flux metric, denoted as Phytotoxic Ozone Dose (PODy) has been adopted by the UNECE Convention on Long-Range Transboundary Air Pollution (LRTAP) to develop doseresponse relationships for the derivation of 'critical levels' for Europe; these are levels below which crop damage would not be expected to occur according to current knowledge (LRTAP Convention, 2017).These 'critical levels' have been used to establish national and regional air quality standards for the formulation of emission reduction policy (Massman et al., 2000;Emberson et al., 2000;Mills et al., 2011).Current flux-response relationships have been developed using an empirical multiplicative g sto model (LRTAP Convention, 2017), a component of the DO 3 SE O 3 deposition model used in European scale modelling (Simpson et al., 2012) to calculate stomatal O 3 flux for crops grown in European filtration/fumigation experiments.This approach allows accumulated stomatal O 3 flux (PODy) to be calculated over a growing season and plotted against relative yield loss for a range of experimental O 3 treatments.A response relationship can then be derived from statistical linear regression of these pooled data points (Pleijel et al., 2022).In Europe, flux-response relationships for wheat are based on data from 4 European countries, encompassing 14 years and 6 cultivars (LRTAP Convention, 2017).
An important criticism and limitation of existing flux-response relationships is that the estimate of g sto is not related to the plant's main physiological requirement for gas exchange, which is the uptake of CO 2 for carbon assimilation by photosynthesis.This creates a disconnect between O 3 stomatal uptake and critical physiological processes such as photosynthesis, respiration, carbon accumulation, and allocation, development, growth, and yield (Ball et al., 1987;Wang et al., 2009).Stomatal conductance models coupled to photosynthesis were developed in the early 1990s (Leuning et.al., 1995) and work on a supply and demand basis whereby stomatal opening regulates the CO 2 availability (supply) and the photosynthetic process in the leaf's chloroplasts determines the plant's need for CO 2 (demand), thereby controlling g sto according to the requirements for photosynthesis.These models are more complex than the empirical multiplicative g sto model since they require an estimate of photosynthesis, which often involves applying a biochemical model to simulate plant physiological processes (Büker et al., 2007;Op De Beeck et al., 2010).However, using a multiplicative model requires more parameters and cannot consider the interaction of different environmental variables at the same time.Using an A net g sto approach would also allow a more mechanistic representation of O 3 effects on growth and yield to be explored (Emberson et al., 2018;Büker et al., 2007).This is important as O 3 is thought to cause damage via both an instantaneous effect on photosynthesis as well as a longer-term effect that induces early onset senescence which may lead to earlier maturity and a shorter time period for grain filling (Ewert and Porter, 2000;Emberson et al., 2018).
In this paper, we develop leaf level A net g sto models suitable for quantifying stomatal O 3 flux.The aims of this paper are (i) to assess the ability of the multiplicative g sto and A net g sto models (an empirical A net g sto model (A net g sto emp) and a mechanistic A net g sto model (A net g sto mech)) to simulate g sto (and A net ), (ii) to assess the ability of A net g sto models to simulate O 3 damage to photosynthesis and leaf senescence, and (iii) to compare the ability of multiplicative g sto and A net g sto models to simulate yield loss and hence derive flux-response relationships.This will be achieved by re-analysis of the European wheat flux-response data used to derive the current UNECE LRTAP Convention flux response relationship (LRTAP Convention, 2017) along with additional data from the UK and Sweden which provide further insight into the effects of O concentrations on leaf physiology and senescence.These three models were not designed to simulate dynamic crop growth or yield but rather to estimate cumulative stomatal O 3 flux for regression against yield to develop flux response relationships.The models can be tested against observed A net , g sto and Chlorophyll Content Index (CCI) data to assess their ability in simulating key aspects of leaf physiology that determine O 3 uptake and damage.

g sto emp model
The g sto emp model is an empirical model that estimates g sto according to environmental modifications to a species-specific maximum stomatal conductance value (g max ) (Jarvis, 1976;Emberson et al., 2000;Pleijel et al., 2007) written as; [ min Where g sto is the flag leaf stomatal conductance (mmol O 3 m − 2 PLA s − where PLA is the projected leaf area) and g max is the species-specific maximum g sto .The parameters leaf f phen , f O3 , f light , f temp , f VPD , and f SWP account for the effect of phenology, O 3 , light, temperature, vapour pressure deficit (VPD), and soil water potential (SWP) on g max .f min is the fractional minimal daylight g sto .These functions have values ranging from 0 to 1. Since wheat grown in the filtration/fumigation studies was always well-watered, we assume that f SWP equals 1.The DO 3 SE algorithms and parameters for these functions are described in equations S1-S5 and Table S1 respectively after Grünhage et al. (2012) and the LRTAP Convention (2017).

A net g sto emp model
The coupled A net g sto emp model provides a consistent estimate of the exchange of CO 2 (driven by supply and demand of CO 2 for photosynthesis and its products) on consideration of water loss controlled by g sto .The A net g sto emp model consists of a combination of two separate models: a) the mechanistic and biochemical photosynthesis model (Farquhar et al., 1980;Harley et al., 1992) that estimates net photosynthesis (A net ), and b) the coupled A net− g sto model of (Leuning, 1995) that estimates g sto .
The A net model assumes that photosynthesis is limited, according to prevailing environmental conditions, by three different mechanisms: i. rubisco activity (A c ); ii. the regeneration of ribulose-1,5-bisphosphate (RuBP) which is limited by the rate of electron transport (A j ) and iii. the rate of transport of photosynthetic products (A p ) (Sharkey et al., 2007).These influences on A net are calculated by determination of the smaller of these theoretical CO 2 assimilation rates less the rate of dark respiration (R d ) (Harley et al., 1992), see equations [2] to [5]. where; Where V cmax25 is the maximum rate of RuBP carboxylation catalysed by the enzyme Rubisco at 25  (Buker et.al., 2007).c s is the CO 2 concentration at the leaf surface and Γ is the CO 2 compensation point, calculated according to Buker et al. (2007).
In the photosynthetic model by Sharkey et al. (2007), the parameters 'a' and 'b' reflect conservative estimates for the electron transport rate during carboxylation and oxygenation, assumed to be 4 and 8 electrons respectively, allowing for the regeneration of RuBP and the formation of NADPH and ATP in the Calvin cycle.A c is modified to include f O3 and leaf f phen to empirically define the effect of leaf age and O 3 induced senescence on g sto (Ewert et al., 1999).This allows V cmax25 to change throughout the growing season.Since O 3 primarily causes a limitation to Rubisco (Ewert et al., 1999), we do not include O 3 damage in estimates of A j and A p .
g sto is calculated from A net using an empirical relationship between g sto , A net and environmental variables following an approach first developed by Ball et al. (1987) and modified by Leuning (1995) as described in equation [6].
Where g min is the minimal daylight g sto value (Leuning,1995).The parameter m describes the species-specific sensitivity to A net and CO 2 concentration at the leaf surface.c s is the CO 2 concentration at the leaf surface and Γ is the CO 2 compensation point calculated according to Buker et al. (2007).The use of the multiplicative g sto models f VPD relationship (Danielsson et al., 2003;Pleijel et al., 2007;LRTAP Convention, 2017) ensures consistency between the g sto emp and A net g sto emp modelling methods used in this study, see equation [7].
where VPDo is the VPD threshold (Leuning et al.,1998) parameterised to reflect a more gradual decrease in g sto with increasing VPD compared to that previously suggested by Leuning's (1995) hyperbolic function (see Fig S1).The A net g sto emp model follows the same method as used in the g sto emp model to calculate the O 3 (i.e. the f O3 function) and phenology (i.e. the leaf f phen function) effect on conducatance.The only structural difference between the A net g sto emp and A net g sto mech model lies in a more mechanistic approach in the latter to model these effects..

A net g sto mech model
The A net g sto mech model simulates the loss of instantaneous photosynthetic activity and the acceleration of leaf senescence using a mechanistic approach to modify the Rubisco-limited rate of photosynthesis (A c ) following the approach of Ewert & Porter (2000) as described in equation [8].
The short-term impact of O 3 on A c is calculated according to the fO 3,s (d) term, the cumulative daylight hour effect of O 3 on Vc max , which allows for an instantaneous effect of O 3 on photosynthesis when stomatal O 3 flux overwhelms detoxification and repair mechanisms (Betzelberger et al., 2012;Feng et al., 2022).fO 3,s (d) is estimated by calculating f O3,s (h) (representing the linear relationship between stomatal O 3 flux (f st )) and a decrease in A c calculated for every hour as described in equation where γ1 and γ2 are both short-term O 3 damage coefficients, with γ1 γ2 representing the O 3 detoxification threshold below which no damage occurs to the photosynthetic system and γ2 determines the effect of f st on A c , see Section 2.2 for the f st calculation, which is estimated for the previous hour.fO 3,s (d) and f O3,s (d − 1) are calculated as described in equation [10].
Where the term f O3,s (d) describes the instantaneous O 3 effect on Vc max25 which is allowed to build over the course of the daylight period (when photosynthetically active radiation (PAR) is greater than 50 W m − 2 ) from an initial value which is determined by the previous days f O3,s (d − 1) value and an allowance for incomplete overnight recovery in Vc max25 which varies with leaf age as described by r O3,s term in equation [11].
Where f LA defines leaf age and is calculated as The long-term impact of O 3 on Vc max25 represented by the f Ls term represents the longer-term accumulation of stomatal O 3 flux (acc fst ) causing degradation to the Rubisco enzyme triggering early and enhanced senescence of mature leaves (Gelang et al., 2000;Osborne et al., 2019).The simulation of f Ls (and f LA used in the short-term O effect) are related to thermal time defined periods over the course of the flag leaf life span defined as a mature (tl,ep) and a senescing (tl,se) stage which together comprise the full flag leaf lifespan (tl, ma), equivalent to leaf f phen in the empirical models.The tl, ep stage defines the period between the start of anthesis and start of senescence (SOS).The tl, se stage simulates the decline in chlorophyll content and depicts the period between SOS and the end of senescence (EOS), see Section 2.4 for the SOS and EOS calculation.TT leaf represents the cumulative thermal time.This value is determined by integrating daily mean temperature over a 24-hour period and accumulating over the course of the growing season.
Equations S5 and S6 give the leaf f phen and tl, ma equations and Fig. S2 describes the relationship between leaf f phen , f LS and f LA .The O effect on f Ls is first simulated by estimating a weighted accumulated fst (fO3, l) modified from Ewert and Porter (2000) by where γ3 determines the reduction in tl, ma as POD y (in µmol m − 2 ) increases and POD y is calculated as described in equation [19].
The SOS is determined by γ4, whilst γ5 determines maturity (or EOS).
Where, tl epO3 is tl ep with an O 3 effect which may bring the onset of senescence earlier, and tl seO3 is tl se with an O 3 effect which may bring maturity earlier.f Ls is estimated by,

Estimation of O 3 uptake (f st ) and PODy
For all models used in this study f st (in nmol O 3 PLA m − 2 s − 1 ) is calculated as a function of O 3 concentration at the leaf boundary layer, g sto and O 3 deposition to the external leaf surface (see equations [17], [18] and [19]) following the LRTAP Convention (2017).
Where [O3] is the O 3 concentration at the upper surface of the quasilaminar boundary layer of the flag leaf (nmol/mol); gsto is leaf stomatal conductance (m/s) as described in Eqn 1 and 6, leaf rb is the quasi laminar leaf boundary layer resistance (s/m), Lm is the cross wind leaf dimension (m), uh is the windspeed at the canopy surface (m/s), leaf rc is leaf surface resistance (s/m), and gext is the external plant cuticle conductance (m/s).Here we assume that the O 3 concentrations measured within the field chambers of the filtration/fumigation experiments represent a reasonable estimate of O 3 at the leaf boundary layer due to the enhanced air circulation.Parameter values are provided in Table S3.This study uses the PODy stomatal flux-based index currently used by the LRTAP Convention (2017) to assess damage to European wheat calculated using a y threshold value of 6 nmol O 3 m − 2 PLA s − 1 according to equation [20] for all three models.
where f sti is the hourly mean O 3 flux in nmol O 3 m − 2 PLA s − 1 (see equation [17]) and n is the number of hours within the accumulation period.y (equivalent to γ1 γ2 ) is equal to 6 (nmol m − 2 PLA s − 1 ) and is subtracted from each hourly averaged f st (nmol O 3 m − 2 PLA s − 1 ) value only when f st > y, during daylight hours (i.e. when PAR > 50 W m − 2 ).The term (3600/10 6 ) converts to hourly fluxes and to mmol O 3 m − 2 PLA.This method estimates POD6 on a per m 2 basis representative of the flag leaf only; it takes no account of the actual LAI of the flag leaf (or other canopy leaves), that might be contributing to carbon assimilate and hence influence O 3 damage.This assumption may warrant further investigation were canopy O 3 uptake considered an important determinant of ozone damage.However, at least for wheat, the importance of the flag leaf in providing carbon assimilate for grain filling likely makes this a reasonable assumption.

Datasets
The g sto models were applied to simulate POD 6 for O 3 filtration/ fumigation experimental datasets conducted since the 1980s in Europe that described wheat yield losses due to different O 3 treatments.These datasets represent 4 countries (Belgium, Sweden, Finland, and United Kingdom) 6 cultivars and 14 years.These are predominantly the same data used to derive the UNECE LRTAP flux-response relationships (LRTAP Convention, 2017) (exceptions being the exclusion of an Italian dataset which used a variety of Durum wheat), and the inclusion of new data from the UK and Sweden which have the benefit of also providing important physiological and chlorophyll content data.A detailed description of these datasets is given in the Table S2.

Parameterisation for the g sto models
The multiplicative g sto model uses the same parameters as described in the LRTAP Convention (2017).Full details are provided in Table S1.
Both the A net g sto emp and A net g sto mech models require parameterisation of V cmax25 , J max25 and m.Parameters, such as g min , representing the minimum stomatal conductance (set to 0.01 µmol CO 2 m − 2 s − 1 ), are sourced from (Ewert and Porter, 2000), while VPD 0 (set at 2.2 kPa and detailed in Section 2.2) are determined empirically.
However, the A net g sto mech model requires additional parameterisation for the O 3 damage module (represented by γ coefficients).By contrast, the A net g sto emp model uses the same fO 3 function as the multiplicative g sto model for estimating O 3 damage and therefore does not need additional calibration.
A systematic literature review was conducted to extract data to define the likely range and initial values (range mean) of V cmax25 , J max25 and m values occurring in wheat across Europe (see section SF); this approach is similar to that used to parameterise the g sto emp model (LRTAP Convention, 2017).V cmax25 and J max25 values were recorded for fully developed flag leaves growing under ambient atmospheric concentrations of O 3 and CO 2 for crops grown in the field/or large pots under a stress-free environment (see Fig. S3).Information describing the bio-geographic region and the prevalence of rainfed or irrigated management were also recorded (Fig. S4).A diagrammatic representation of the systematic literature review is provided in Fig. S5.
The parameterisation of m needs to be considered in relation to VPD 0 since the slope of the relationship m found when plotting A net against g sto represents a compromise between the cost and benefit of g sto relative to CO 2 uptake for photosynthesis vs water loss affecting intrinsic water use efficiency (Medlyn et al., 2011).Here we follow the approach of Medlyn et al. (2011) and calibrate m to ensure that the modelled maximum A net against g sto aligns with the maximum observed A net against g sto values.
The parameters γ3, γ4, and γ5 are only used in the A net g sto mech damage module to simulate the rate of senescence.They were calibrated to ensure that the start (SOS) and end (EOS) of the senescence period matched observed senescence timings.These observations were derived from data describing the Chlorophyll Content Index (CCI) using the 'break point' analysis method (Mariën et al., 2019).This method determines the change in the seasonal pattern of CCI (and hence senescence) as a function of day of the year through piecewise linear regressions.The first segment of the regression (i.e.leaf expansion to mid-anthesis) was constrained to zero since it is assumed the leaf does not undergo senescence during this period.The slope of the second segment (from mid-anthesis to harvest) was allowed to be greater than zero on the assumption that senescence of the flag leaf will only occur after mid-anthesis.The slope with the lowest RMSE, indicating the smallest deviation between the measured CCI data points and the values estimated by the piecewise linear regression model, was assumed as the breakpoint for the SOS.Furthermore, a polynomial regression line, which delineates the period of senescence, was employed to determine EOS.The SOS and EOS of the flag leaf determined from break-point P. Pande et al. analysis of the UK (2015) and Swedish (1997 and 1999) datasets are given in the Table S4.
Details of the initial values and associated ranges for calibration of all A net g sto parameters are provided in Table 1.

Calibration of the A net g sto emp and A net g sto mech models
The parameters for the g sto emp model were taken directly from LRTAP Convention (2017) and as such further calibration adjustments were not performed in this study.
The A net g sto mech and A net g sto emp model calibration for European conditions is performed in steps (as outlined below) using g sto , A net and CCI data from various sub-sets of the fumigation/filtration dataset.Fig. 1 presents a schematic diagram of the calibration process used for the A net g sto models.
In the first step, initial values for V cmax25 , J max25 and m are selected that give a maximum g sto value of between 500 and 600 mmol O 3 m − 2 PLA s − 1 and a maximum A net value of between 30 and 35 µmol CO 2 m − 2 s − 1 .These values are consistent with the experimental dataset for Bangor as well as published studies that provide values for these parameters across Europe (Uddling and Pleijel, 2006;Sharma et al., 2015).This step only uses the low O 3 treatment data from Bangor (n = 14, see section SH) to ensure leaf physiology is unaffected by O 3.
In the second step, which is only performed for the A net g sto mech model, the focus is on establishing initial values for O 3 damage parameters (γ1 to γ5) using datasets from both low (n = 11) and very high (n = 10) O 3 treatments from Bangor (see section SH).The O 3 coefficients γ1 and γ2 were set to give a detoxification threshold of 6 nmol O 3 m − 2 s − 1 , while γ3, γ4, γ5 were calibrated based on the observed SOS and EOS data, identified using the breakpoint method discussed in Section 6. O 3 damage parameters for the A net g sto emp model are used as provided in the LRTAP Convention (2017) based on the f O3 function (and so consistent with the methods used in the g sto emp model).
Moving to the third step, model calibration uses all O 3 treatment data, segmenting these data into training and test sets as detailed in the Table S5.This uses a bootstrapping resampling technique (Hesterberg, 2011), using R software 4.2.3, to create bootstrap samples (n = 5) that randomly select a dataset with replacement i.e., in a sample, there can be duplicates of the same dataset (Table S5).Such an approach ensures that the initial parameters from steps one and two, along with their defined ranges drawn from both these steps and existing literature, are robustly tested across diverse data combinations from the fumigation/filtration experiments.
The calibration process then proceeds with these training samples (n = 5), aiming to calibrate the model to find the best parameters for V cmax25 , J max25 and m, and O 3 damage parameters (γ3 to γ5, only for the A net g sto mech model).This calibration employs a computational genetic algorithm (Wang, 1997), an optimisation technique, with gradient descent to find the best parameters.The process requires an initial value and a range, and 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 with the highest R 2 and lowest RMSE value.Finally, the calibration outcomes from each training sample are aggregated, using weighted averages following Eq.S7, to establish the final  set of parameters.These parameters are then used to run the models to estimate POD 6 and hence construct the flux-response relationships (Fig. S7), ensuring the model's applicability and accuracy.The model's efficacy is then tested using test datasets (n = 5), which apply these final parameters.The performance metrics for these tests, specifically the R 2 and RMSE values for the flux-response relationships, give an indication of the model's reliability and precision across different datasets.

Leaf physiology
Leaf physiology data (g sto and A net ) from the UK were used to assess the ability of the different models to simulate key physiological variables necessary to estimate POD y under both low background and peak O 3 treatments over the course of the growing season.
Fig. 2a and b show a scatter plot of model simulations of hourly mean g sto values plotted against observed values for the 2015 and 2016 background and peak O 3 treatments for Mulika and Skyfall wheat varieties.All g sto models performed similarly under the background O 3 treatments with R 2 values of between 0.33 and 0.43 and RMSE values between 111 and 137 mmol O 3 m − 2 s − 1 with the A net g sto mech model performing the best.All g sto models performed less well under the peak O 3 treatment with the R 2 range between 0.07 and 0.33, with the A net g sto mech model performing the best; all models have similar RMSE values.For the peak O 3 treatment, the A net g sto mech model tends to overestimate g sto whilst the other two models tend to underestimate g sto in relation to the 1:1 line.Similar results were found for A net with values simulated reasonably well under background O 3 treatments by both the A net g sto emp and A net g sto mech models with R 2 values of between 0.8 and 0.83 (see Fig. S8a).Both the models tend to underestimate maximum values of A net by ~10 µmol CO 2 m − 2 PLA s − 1 .
All models were able to simulate the mean diurnal (see Fig. S9) and mean daily maximum (see Fig. 2c) g sto values equally well for the background O 3 treatment.For the peak O 3 treatments, the A net g sto mech model tended to overestimate mean diurnal g sto by about 50 mmol O 3 m − 2 PLA s − 1 whilst the other two models tended to underestimate g sto by the same margin.Similarly, models were able to simulate the mean diurnal (see Fig. S10) and mean daily maximum A net values (see Fig. S8c) equally well for the background O 3 treatment.As for g sto , all models struggled to predict A net under the peak O 3 treatments with a tendency to overestimate A net in relation to the 1:1 line but to underestimate maximum A net values.A net was comparatively better predicted by the A net g sto mech model with R 2 values of 0.42 compared to 0.31 for A net g sto emp model.
Fig. 2c shows that the A net g sto mech model performs better under peak O 3 concentrations over the full length of the flag leaf lifespan, thus simulating the effect of senescence on g sto reasonably well.By contrast the g sto emp and A net g sto emp models simulated an overly sensitive senescence response of g sto to O 3 compared to the observations.Similar to the g sto results, the A net g sto models overestimated the decline in A net at the end of the growing season compared to the observations (see Fig. S8b).However, the A net g sto mech model gave a closer fit to the observations than the A net g sto emp model.It is also worth noting that the A net g sto mech model simulates higher g sto and A net under the peak O 3 treatment than the low O 3 treatment for the UK.This is because the O 3 effect is most strongly determined by its longer-term impact on senescence than its instantaneous impact on photosynthesis, the former only taking effect once O 3 has brought forward the SOS which occurs only towards the end of the growing season where there are far fewer observed data for comparison.

Leaf senescence
The CCI data available from the UK (cv Mulika) and Swedish (cv Dragon) filtration/fumigation datasets were used with the break point method to estimate the SOS and EOS.Results in Fig. 3 show that the higher O 3 treatment (low background vs very high peaks for the UK data) brought forwards the SOS by 7 days and EOS by 12 days.Similar results are found for Sweden by comparing the CF vs NF++ experiment with SOS and EOS being brought forwards by 6 days and 12 days respectively (see Fig. S11).
The data provided in Table 2 can be used to assess the ability of the A net g sto models to simulate senescence under the different datasets and O 3 treatments used in this study.Table 2 summaries information for the extreme O 3 treatments (i.e.comparing lowest with highest).The difference in O 3 treatment causing senescence effects is indicated by the POD 6 values for the flag leaf lifespan.Table 2 shows that the A net g sto emp model predicts SOS to occur earlier with a range of 20 days difference compared to the observations, and EOS to generally occur later with a range of 18 days difference compared to the observations.By comparison the A net g sto mech model simulates SOS closer to the actual date with a range of 8 days earlier to 3 days later and EOS with a range of 8 days earlier to 3 days later compared to the observations.The POD 6 values for the high O 3 treatments are consistently higher for the A net g sto mech model suggesting that the model is paramterised to be less sensitive to cumulative stomatal O 3 uptake than the A net g sto emp model.Overall, the mechanistic approach used by the A net g sto mech model simulated SOS and EOS more closely to the observations.However, care should be made in interpreting these results since the CCI data used to define the actual SOS and EOS are limited in number, leading to some uncertainty in the actual timings of senescence, especially close to anthesis.It should also be noted that the A net g sto models are calibrated against all the CCI data held in the datasets and so there will be some discrepancy when comparing simulations against individual datasets and O 3 treatments.

Flux-response relationships
Each of the three g sto models were used to develop flux-response relationships based on POD 6 using the O 3 filtration/fumigation data (Fig. 4).The robustness of the flux-response relationship can be determined by the strength of the linear regression (i.e., R 2 value).The A net g sto mech model (R 2 = 0.74) performed better than the g sto emp model (R 2 = 0.68) in deriving flux-response relationships.The A net g sto emp model performed slightly less well (R 2 = 0.66).The slope of the relationships differ by − 0.0412, − 0.0342 and − 0.0325 for g sto emp, A net g sto emp and A net g sto mech respectively.This is because the A net g sto mech model simulates higher g sto values under elevated O 3 and during senescence which will increase the POD y values.This demonstrates the importance of consistency in using the same g sto method to estimate POD y as is used to derive the flux-response relationship for yield loss estimates.Were 'critical levels' to be derived from these relationships using the methods described in the LRTAP Convention (2017) (i.e. a 5% reduction in grain yield based on the slope of the relationship) values of 1.69, 1.19 and 1.75 mmol O 3 m − 2 would be found for g sto emp, A net g sto emp and A net g sto mech models respectively (also shown as dotted lines in the respective plots in Fig. 4).The range of these values reflects the high g sto values modelled using the A net g sto mech model.It is useful to note that the dose-response relationships developed in this study are an improvement to those presented in the LRTAP Convention (2017) Mapping Manual (albeit with slightly different data compliments).For comparison, we also show the dose-response relationships developed by applying these three models but only with those datasets used in the LRTAP Convention (2017) Mapping Manual (see Fig. S12).

Discussion
We found that the process-based A net g sto mech model can derive robust flux-based dose-response relationships (with an R 2 value of 0.74), this performance is marginally improved to that of empirical-based  In each plot, the red solid line represents the regression line, showing the relationship between the modelled and observed values.The black dashed line represents the 1:1 line, the coefficient of determination (R 2 ) and root mean square error (RMSE) is provided for the regression; and b) Average daily maximum g sto values simulated over the flag leaf lifespan by each of the three g sto models and observed daily maximum g sto data.Standard error bars for the observed data are given by black lines extending from the red observed points, providing a visual representation of uncertainty in the measurements.models (g sto emp and A net g sto emp).However, there was little difference (ranging from 1.69 to 1.75 mmol O 3 m − 2 ) in the 'critical levels' derived using each method.This suggests A net g sto models can be reliably used in the derivation of dose-response relationships and 'critical levels' for regional scale risk assessments and that the slope of the dose-response relationship was robust from the point of view of the method to model POD 6 .However, even though the variability in slope and 'critical level' values are relatively small, these differences highlight the importance of consistency in application, i.e. that the same g sto algorithm be used to derive the flux (POD y )-response relationship used in the risk assessment.Our study also found that the A net g sto mech model was better able to simulate the diurnal and seasonal variation in observations of both A net and g sto found under low vs high O 3 treatments in the Bangor experiment.This model attribute is particularly advantageous in estimating POD y given that O 3 concentration profiles can vary substantially across the global wheat growing regions, with some experiencing more chronic O 3 concentrations (e.g., in Europe (Karlsson et al., 2017) while others will experience more extreme, episodic concentrations (e.g., in Asia (Lei, Wuebbles and Liang, 2012)).The results suggest that the A net g sto mech is better able to simulate stomatal O 3 uptake under conditions of higher O concentration.Since the slope of the resulting dose-response relationship does not change, this suggests that the sensitivity of wheat to O uptake remains consistent but that the model is better able to simulate what actual uptake occurs.This finding would warrant further

Table 2
Comparison of the difference in days between Start (SOS) and End (EOS) of senescence by site, year and O 3 treatment (described by average 24-hour mean O concentrations in ppb).The "SOS bias" and "EOS bias" columns indicate the deviation in days at SOS and EOS, respectively, from applying the A net g sto models as compared to the observed data.Positive values denote a delay, while negative values signify an advancement in the modelled timing of senescence relative to the observations.Also shown are the POD y values at SOS and EOS.investigation as new datasets become available.
There are three important aspects to accurate A net and g sto estimates, firstly the parameterisation of the leaf level A net model which is dependent upon V cmax25 , J max25 , m and VPD 0 .Secondly, the instantaneous effect of O 3 on A net in relation to its parametrisation and effectiveness in causing O 3 damage.Thirdly, the parameterisation of the module describing O 3 induced leaf senescence, the latter is especially important to estimate A net and g sto toward the end of the growing season, in wheat this coincides with the grain-filling period and is therefore important in determining yield (Neghliz et al., 2016).
Parametrised values for V cmax25 and J max25 of 88 and 173 µmol CO 2 m − 2 s − 1 respectively in this study compare reasonably well to the values of 62-75 and 150-195 µmol CO 2 m − 2 s − 1 used for LINTULLC2 (Feng et al., 2022) and AFRCWHEAT (Van Oijen and Ewert, 1999) crop models which incorporate O 3 damage modules for similar European wheat applications.We found limited evidence for variation in V cmax25 and J max25 with biogeographical region with V cmax25 varying between 55 and 180, 53 and 185 and 90 and 120 µmols CO 2 m − 2 s − 1 for Atlantic, continental and Mediterranean biogeographic regions respectively; no statistical difference by region was found.This contrasts with the g max value of the g sto emp model that has lower values for Mediterranean wheat cultivars (by 70 mmol O 3 m − 2 s − 1 , LRTAP Convention ( 2017)).This study only used experimental data from Atlantic, Boreal or Continental regions.Were Mediterranean data to have been included, the V cmax25 and J max25 values may have warranted further investigation to establish whether a different Vc max25 might be justified, especially since only 11 datapoints were retrieved for this region in our literature search (see Fig. S3).An indication of this can be provided through comparison with the modelling study presented by Nguyen et al. (2024) Nguyen et al. (2024) which used three crop models (including DO 3 SE-Crop, an extension of the A net g sto type of model described here to estimate carbon allocation, growth, and yield).Here the DO 3 SE-Crop model was parameterised for a Mediterranean variety of spring wheat (Califa sur) with values of key photosynthetic model parameters being 102 µmol CO 2 m − 2 s − 1 for V cmax25 , 194 µmol CO 2 m − 2 s − 1 for J max25 , 8.57 for m and 2.2 kPa for D 0 .These Mediterranean values for V cmax25 and J max25 are both somewhat lower (by ~14 and 20 µmol CO 2 m − 2 s − 1 respectively) than those used in this study and hence would suggest that the A net g sto model would benefit from a Mediterranean parameterisation similar to the regional parameterisations used in the g sto emp model.
One other important consideration in relation to geographical region is the effect of soil moisture on g sto since the Mediterranean region is likely to experience longer and more extreme periods of drought stress that will reduce stomatal O 3 uptake (Fagnano et al., 2009).This is particularly important for wheat since this tends to be a rainfed crop in Europe.A variety of methods have been developed to simulate the effect of soil water status (described variously as soil water potential (referred to in this study as f swp ), soil water content or plant available water (Büker et al., 2007) on g sto .These methods can be used in either the g sto emp or A net g sto emp type models (the latter by including the f swp function as a multiplier in the A c formulation (see Eq. ( 3)).We were unable to test the effectiveness of this aspect of the modelling since the datasets used in this analysis all represented well-watered conditions.However, this would be an important aspect to investigate further, especially in relation to model application, to ensure g sto emp and A net g sto emp models respond similarly (in terms of magnitude of changes to stomatal O 3 flux) to the inclusion of these soil water status parameters.
The ratio between V cmax25 and J max25 was found to vary between 0.2 and 0.8 (Fig. S3) and was calibrated to a value of 0.51 for this dataset.This is consistent with a study by Wullschleger (1993) who found a ratio of 0.38-0.55 for wheat even as growth and temperature varied.However, other research found that the ratio may range from 1 to 3 (Camino et al., 2019;Day,Station and Al, 1982) which may be attributed to J max25 being more reliant on light than V cmax25 causing the ratio to decrease when light intensity decreases (Dai et al., 2004).The value of 7.87 for m used in this study is also within the range of 5 and 15 found for many different cultivars of wheat (Kosugi et al., 2003;Collatz et al., 1991;Baldocchi and Meyers, 1998;Miner,Bauerle and Baldocchi, 2017).The VPD 0 value is markedly different (2.2 kPa) from that of Luening et al. (1995) and means that A net can be maintained under high values of VPD, this is consistent with the f VPD relationship and observational data (Danielsson et al., 2003).
The validity of the A net g sto mech model also depends on appropriate formulation and parameterisation of the key O 3 damage mechanisms.These damage mechanisms are assumed to have both an instantaneous (f O3,s (d)) effect of O 3 on photosynthesis and a longer-term effect (f LS ) of accumulated O 3 uptake promoting earlier senescence.The instantaneous effect reduces carboxylation via a reduction in rubisco activity which may in turn lead to a reduction in carbon assimilation when Rubisco activity ( A c ) is limiting net photosynthesis.This reduction in Rrubisco activity is assumed to repair overnight but with repair effectiveness decreasing as the leaf ages.According to Farage et al. (1991), the instantaneous impact of O 3 was only seen with a significant reduction in carboxylation efficiency (>50 %) causing a reduction in carbon assimilation.This could happen when crops are exposed to elevated O concentrations for long periods or if repeated high O 3 exposures were to take place causing the crop to lose its ability to recover (Feng et al., 2022).By contrast, the length of the leaf senescence period is essential for determining the crop development cycle (Ding et al., 2023).The onset of leaf senescence causes a substantial decrease in carbon assimilation (A net ), primarily attributed to changes in chloroplast structure and function, and hence the chlorophyll content in the flag leaf (Ding et al., 2023;Gelang et al., 2000;Ojanperä et al.,1998), and contributes to the reduction in dry ear weight, which directly affects yield loss (Gelang et al., 2000).The CCI has been shown to be a good predictor of the onset of senescence (Mariën et al., 2019;Osborne et al., 2019).It can also be used as a proxy for V cmax25 (Croft et al., 2017), which is our modelling approach since we assume SOS will coincide with a reduction in V cmax25 and consequently A c (see Eq. ( 8)).We find that the A net g sto mech model can simulate SOS and EOS for the elevated O 3 treatments in the UK and Sweden data better than the empirical models.For the UK, the flag leaf starts to senesce 6 days earlier in high (VHP) compared to low (LB) O 3 treatment, for Sweden 7 days earlier in high (NF+++) compared to carbon-filtered (CF) treatments.The number of days by which high O levels can bring forward the start of senescence is corroborated by other published studies (Pleijel et al., 1997;Grandjean and Fuhrer, 1989;Gelang et al., 2000) which found the flag leaf could senesce up to days earlier in the very high O 3 compared to the carbon filtered treatments.O 3 was also found to cause differences in the maturity (EOS)of the flag leaf; Shi et al. (2009) reported that maturity (EOS) occured days earlier in elevated O 3 (50 % higher than ambient) compared to ambient O 3 treatments.Similar results were found in this study, with the flag leaf modelled to reach maturity (EOS)12 days earlier in VHP compared to LB treatments.Although our results seem consistent, they are based on a limited number of CCI data points (11 and 13 for each treatment for the UK and Sweden respectively) which are only captured from mid-anthesis to 10 days before maturity.Additional CCI data spread more evenly over the crucial crop growth period would improve our understanding of how O 3 affects senescence.
Parameters for the A net g sto models were found using an automated calibration method, the genetic algorithm optimisation technique since this approach is considered superior in performance to more traditional techniques (Kuo et al., 2000;Dai et al., 2009;Vazquez-Cruz et al., 2014).The genetic algorithm method was also chosen since it works with a range of parameter searches from a population of points and employs probabilistic transition rules, i.e., uses random sets of parameters instead of using fixed sets, which makes the optimisation process more robust (Kuo,Merkley and Liu, 2000).This study demonstrated the effectiveness of this approach with the five training samples that are used to form dose-response relationships giving RMSE ranges from 0.99 to 4.5 × 10 − 5 mmol m − 2 for the A net g sto mech model (data not shown).Such a good performance suggests that the parametrisation derived can give robust values for the A net g sto models for use in other European O 3 risk assessment applications.
The calibration approach to parameterise the A net g sto models is different to that used to parameterise the g sto emp model which identifies g max and f min values (as average maximum and minimum values respectively) and the f functions using a boundary line analysis method (LRTAP Convention, 2017).Since the A net models are effectively calibrated to the output of a sub-set of all datasets it can be argued that this may improve the ability of this model type compared to the g sto emp model.It is also important to note that the A net models calibration included the UK Bangor dataset (and hence additional information on the onset and rate of senescence) as compared to the parameterisation of the g sto emp model, these data would have been useful to test and inform the existing g sto emp f O3 function.Ideally, all models would be calibrated using the same data and methods, which would mean that the g sto emp model would be calibrated using the genetic algorithm method and with the inclusion of the UK data describing senescence.Although such work was outside the scope of the current study it would be useful to consider in future modelling studies.As such unequivocal claims that A net g sto models are better than g sto emp models need to be made with caution.

Conclusion
Overall, we find that the A net g sto mech model can be used to derive robust flux-response relationships when incorporating both short-and long-term O 3 damage processes.The A net g sto mech model also has the added benefit of achieving reasonable estimates of g sto under variable O 3 concentrations and has a direct link to carbon assimilation.This study's establishment of an A net g sto mech flux-response relationship could be used to calibrate or constrain models that use the A net g sto approach (e.g.photosynthesis based crop models, land surface exchange models, biogeochemical cycling models and earth system models) thus supporting a move towards more process-based assessments of O 3 damage and yield loss.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Schematic diagram of the calibration process used for the A net g sto models.This describes the number of filtration/fumigation datasets used for both the training and testing of model performance in relation to the automated calibration of various parameters dependent upon the construct of the A net g sto emp and A net g sto mech models.'n' and 'z' refer to the number of datasets and parameters used, respectively.

P
.Pande et al.

Fig. 2 .
Fig. 2. Plots for background and peak O 3 treatments for Mulika and Skyfall wheat cultivars, fumigated in Bangor over the 2015 and 2016 growing seasons showing a) Observed against modelled g sto values estimated using the three different g sto models.In each plot, the red solid line represents the regression line, showing the relationship between the modelled and observed values.The black dashed line represents the 1:1 line, the coefficient of determination (R 2 ) and root mean square error (RMSE) is provided for the regression; and b) Average daily maximum g sto values simulated over the flag leaf lifespan by each of the three g sto models and observed daily maximum g sto data.Standard error bars for the observed data are given by black lines extending from the red observed points, providing a visual representation of uncertainty in the measurements.

Fig. 3 .
Fig. 3. Leaf senescence profiles of O 3 induced leaf senescence for the Mulika wheat cultivar for the low background (LB) and very high peak (VHP) O 3 treatments in the UK dataset.The timing of the SOS and EOS (vertical dotted black lines) determined by applying the break point method to the CCI data (red circle with standard error bars) are shown in relation to estimates made by the A net g sto emp model (which uses leaf f phen and f O3 functions to simulate senescence and the A net g sto mech model (which uses f LS ) to simulate senescence.

Fig. 4 .
Fig. 4. Flux-response relationships for relative wheat grain yield derived using the three g sto models to simulate the POD 6 metric.The plots replicate the LRTAP Convention (2017) dose-response relationships with the exception of exclusion of an Italian for Durum wheat, and inclusion of UK and an additional Swedish dataset.The 95% confidence intervals are indicated by the dotted lines around the best-fit line.The vertical dashed line indicates the 'critical levels' determined by each model.Each figure includes the coefficient of determination (R 2 value) and a dose-response relationship equation.

P
.Pande et al.
• C (leaf temperature), C i and O i are the intercellular CO 2 and O 2 concentrations respectively; K c and K o are the

Table 1
A detailed overview of the parameters, ranges, and optimised values after calibration of the A net g sto models.
*γ parameters only used for A net g sto + O 3 mech .