Modeling the effects of tropospheric ozone on the growth and yield of global staple crops with DSSAT v4.8.0

Elevated surface ozone (O₃) concentrations can negatively impact growth and development of crop production by reducing photosynthesis and accelerating leaf senescence. Under unabated climate change, future global O₃ concentrations are expected to increase in many regions, adding additional challenges to global agricultural production. Presently, few global process-based crop models consider the effects of O₃ stress on crop growth. Here, we incorporated the effects of O₃ stress on photosynthesis and leaf senescence into the Decision Support System for Agrotechnology Transfer (DSSAT) crop models for maize, rice, soybean, and wheat. The advanced models reproduced the reported yield declines from observed O₃-dose field experiments and O₃ exposure responses reported in the literature (O₃ relative yield loss RMSE <10 % across all calibrated models). Simulated crop yields decreased as daily O₃ concentrations increased above 25 ppb, with average yield losses of 0.16 % to 0.82 % (maize), 0.05 % to 0.63 % (rice), 0.36 % to 0.96 % (soybean), and 0.26 % to 1.23 % (wheat) per ppb O₃ increase, depending on the cultivar O₃ sensitivity. Increased water deficit stress and elevated CO₂ lessen the negative impact of elevated O₃ on crop yield, but potential yield gains from CO₂ concentration increases may be counteracted by higher O₃ concentrations in the future, a potentially important constraint to global change projections for the latest process-based crop models. The improved DSSAT models with O₃ representation simulate the effects of O₃ stress on crop growth and yield in interaction with other growth factors and can be run in the parallel DSSAT global gridded modeling framework for future studies on O₃ impacts under climate change and air pollution scenarios across agroecosystems globally.

1 Introduction

Surface or ground-level, ozone (O₃) is a major air pollutant that causes adverse impacts on agricultural productivity worldwide (Mills et al., 2018a; Emberson et al., 2018; Tai et al., 2021). O₃ is formed through photochemical reactions between incoming solar radiation and primary pollutants such as nitrogen oxides (NO_x= NO + NO₂), volatile organic compounds (VOCs), carbon monoxide (CO), or methane (CH₄) across all areas of the globe (Cooper et al., 2014; Simpson et al., 2014). Global O₃ concentrations have increased 2 %–7 % per decade in northern mid-latitude regions and 2 %–12 % per decade in tropical regions since the mid-1990s (IPCC, 2021; Arias et al., 2021). Future O₃ concentrations are projected to continue increasing if O₃ precursor emissions are not mitigated, i.e., following the shared socio-economic pathways where regional rivalry leads to doubling of CO₂ emissions by 2100 (SSP3-7.0) or where fossil-fuel-enabled growth leads to doubling of CO₂ emissions by 2050 (SSP5-8.5) (IPCC, 2021; Arias et al., 2021; Szopa et al., 2021; Griffiths et al., 2021).

Crops exposed to elevated levels of O₃ concentrations can experience reduced photosynthesis, accelerated senescence, foliar chlorosis, and even necrosis from increased cumulative oxidative stress (Ainsworth, 2017). These negative effects lead to decreased productivity resulting in global yield losses between 2 %–16 % for the four main staple crops, maize, rice, soybean, and wheat (Ainsworth, 2017; Schiferl and Heald, 2018; Emberson, 2020), with global annual economic damages of approximately USD 34 billion (Sampedro et al., 2020; Feng et al., 2022). Climate change may exacerbate the negative effects from elevated O₃ concentrations because O₃ concentrations are highest in summer months, and the projected higher temperatures with more frequent heat waves may lead to a longer period of more active photochemical reactions (Zhang and Wang, 2016; Hou and Wu, 2016; Szopa et al., 2021). Elevated concentrations of atmospheric CO₂ and increased periods of water deficit stress cause stomatal closure that can reduce crop O₃ uptake (Khan and Soja, 2003; Biswas et al., 2013), but in turn potential yield gains associated with the CO₂ fertilization effect (Toreti et al., 2020; Jagermeyr et al., 2021) may be constrained by elevated O₃. Therefore, it is important to evaluate net O₃ effects for crop growth and consider the effects of O₃ in global agricultural assessments examining future scenarios.

Process-based crop simulation models have been used to evaluate the impacts of O₃ on crop yields (Guarin et al., 2019; Tai et al., 2021), but most global gridded process-based crop models are still unable to respond to O₃ stress. Recently, the global Lund–Potsdam–Jena managed Land (LPJmL) and Joint UK Land Environment Simulator (JULES) models were modified to include the effects of O₃ stress on soybean and wheat growth (Schauberger et al., 2019; Leung et al., 2020). Additionally, the Agricultural Model Intercomparison and Improvement Project (AgMIP; Rosenzweig et al., 2013) ozone team has recently developed protocols for incorporating O₃ stress into a wider body of crop models aiming to establish the first multi-model assessment of ozone impacts on agriculture at a global level (Emberson et al., 2018).

The aim of this study is to incorporate the effects of O₃ concentrations into the stress response functions of the maize, rice, soybean, and wheat models within the established Decision Support System for Agrotechnology Transfer (DSSAT) v4.8.0 modeling platform (Jones et al., 2003; Hoogenboom et al., 2019, 2021), and consequently the parallel DSSAT (pDSSAT) v4.8.0 global gridded modeling platform that is used to run DSSAT in a global setup (Elliott et al., 2014), to simulate O₃ effects on global crop development and yield for the four major staple crops. The observational data from the free-air CO₂ enrichment (FACE) field experiments conducted in Champaign, Illinois, USA (Choquette et al., 2020; Betzelberger et al., 2012), and well-known O₃ exposure relationships reported in the literature are used to develop and calibrate the model O₃ response functions. Additionally, the observed interactions between O₃, CO₂, and water deficit stress are examined via sensitivity analyses conducted with the modified models.

2 Materials and methods

2.1 Description of crop models

The crop models within the pDSSAT parallel modeling environment are based on the existing crop models within the widely used DSSAT crop modeling platform (Jones et al., 2003; Hoogenboom et al., 2019, 2021) combined with the Center for Robust Decision-making on Climate and Energy Policy (RDCEP) parallel System for Integrating Impact Models and Sectors (pSIMS) framework (Elliott et al., 2014) to allow for global gridded process-based crop modeling on high-performance computational systems. The O₃ stress routines presented here are also applied in the standard DSSAT crop models and can be used for field-level simulations and point-based testing in addition to the global-level modeling applications.

The four DSSAT crop models used in this study are the Crop Environment Resource Synthesis (CERES) -Maize, CERES-Rice, Crop Growth Simulation (CROPGRO) -Soybean, and Nitrogen Wheat (NWheat) models that have been used in previous AgMIP crop model intercomparisons (Bassu et al., 2014; Li et al., 2015; Asseng et al., 2015; Kothari et al., 2022). The CERES-Maize and CERES-Rice models were previously used to estimate global ozone crop losses (Schiferl and Heald, 2018); however, their approach was based on the multiplication of the simulated global base production by the relative yield–O₃ response functions to determine a response proxy. The approach used in this present study integrates daily process-based stress calculations to simulate daily crop growth and stress dynamics. Thus, the models are more applicable to a much broader range of scenarios given that they can combine daily stress interactions and can be used to scale across agroecosystems in a more robust way.

2.2 O₃ incorporation into the crop models

The incorporation of O₃ effects into the DSSAT crop models followed the same methodology as the O₃ incorporation into the DSSAT-NWheat crop model (Guarin et al., 2019), which was based on the incorporation of previous abiotic stress routines (Asseng et al., 2004). O₃ response was added to the models via the inclusion of daily photosynthesis reduction and leaf senescence acceleration functions. Additionally, the interaction between O₃ and water deficit stress and/or atmospheric CO₂ concentrations was incorporated into the models since these combined interactions can mitigate impacts from O₃ on crop production and vice versa. For example, water deficit stress that induces stomatal closure in turn limits O₃ stress because of reduced aerosol uptake (Khan and Soja, 2003; Biswas et al., 2013).

2.2.1 CERES-Maize and CERES-Rice models

The effects of O₃ were incorporated into the CERES-Maize and CERES-Rice models using similar methodology since these two models share similar code. O₃ was added into the models using a photosynthesis reduction stress factor (FO₃) following Eq. (1):

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\text{FO}}_{\mathrm{3}}& =\max \left(\mathrm{0.0},-\left({\displaystyle \frac{{\text{FOZ}}_{\mathrm{1}}}{\mathrm{100}}}\right)\cdot {\text{OZON}}_{\mathrm{7}}\right.\\ & \left.+\left(\mathrm{1.0}+\left({\displaystyle \frac{{\text{FOZ}}_{\mathrm{1}}}{\mathrm{100}}}\right)\cdot \mathrm{25.0}\right)\right),\end{array}\end{array}$$

where OZON₇ is the daily mean 7 h (M7, 09:00–15:59) O₃ concentration (ppb) and FOZ₁ is the O₃ stress parameter for photosynthesis calibrated for different O₃ sensitivities of cultivars divided by a decimal correction factor of 100. The decimal correction factor ensures that the FOZ₁ parameter value ranges between 0.0 and 1.0 in the model ecotype parameter file for comprehensible user input. A minimum M7 O₃ threshold of 25 ppb was set as the reference value based on pre-industrial O₃ concentrations and the United States National Crop Loss Assessment Network (NCLAN) studies indicating that O₃ damage within crops occurs above this threshold (Heck et al., 1984; Lesser et al., 1990; Feng and Kobayashi, 2009). When the daily M7 O₃ concentration exceeds this threshold, photosynthesis is reduced by a factor between 0.0 and 1.0 (Eq. 1), and leaf senescence is accelerated by a factor between 0.0 and 1.0 (Eq. 5). The M7 O₃ metric was chosen as the model input because it is the most readily available metric in the literature, and conversion functions exist to convert between M7 and AOT40, daily mean 12 h (M12), or daily mean 24 h (M24) O₃ metrics (Osborne et al., 2016).

Equation (1) does not include the interaction of O₃ stress with water deficit stress or elevated atmospheric CO₂. To consider these combined interactions on crop growth (PRFO₃), FO₃ was modified using Eq. (2):

$$\begin{array}{}\text{(2)}& {\text{PRFO}}_{\mathrm{3}}=\min \left(\mathrm{1.0},\left({\displaystyle \frac{{\text{FO}}_{\mathrm{3}}\cdot {\text{PCO}}_{\mathrm{2}}}{\text{SWFAC}}}\right)\right),\end{array}$$

where PCO₂ is the atmospheric CO₂ effect on potential daily dry matter production and SWFAC is the water stress factor on photosynthesis (Jones and Kiniry, 1986; Ritchie et al., 1987; Jones et al., 2003). Since PCO₂ is always greater than 1, multiplying by the CO₂ effect mitigates the reduction caused by FO₃. Because SWFAC is a reduction factor between 0 and 1, dividing by this factor decreases the reduction from FO₃ under increased water deficit stress conditions.

The simulated daily biomass production (CARBO, g plant⁻¹ d⁻¹) within the models was calculated based on the existing photosynthesis stress factors with the addition of PRFO₃ using Eq. (3) for maize and Eq. (4) for rice:

$$\begin{array}{}\text{(3)}& {\displaystyle }\begin{array}{rl}{\mathrm{CARBO}}_{\mathrm{maize}}& =\text{PCARB}\\ & \cdot min\left(\text{PRFT},\text{SWFAC},\text{NSTRES},\right.\\ & \left.{\text{PSTRES}}_{\mathrm{1}},\text{KSTRES},{\text{PRFO}}_{\mathrm{3}}\right)\cdot \text{SLPF},\end{array}\text{(4)}& {\displaystyle }\begin{array}{rl}{\mathrm{CARBO}}_{\mathrm{rice}}& =\text{PCARB}\\ & \cdot \text{min}\left(\text{PRFT},\text{SWFAC},\text{NSTRES},\text{TSHOCK},\right.\\ & \left.{\text{PSTRES}}_{\mathrm{1}},{\text{KSTRES,PRFO}}_{\mathrm{3}}\right)\cdot \text{SLPF},\end{array}\end{array}$$

where PCARB is daily potential dry matter production of the crop accounting for light interception, radiation use efficiency, and the CO₂ effect on photosynthesis (g plant⁻¹); PRFT, SWFAC, NSTRES, TSHOCK (CERES-Rice only), PSTRES₁, KSTRES, and PRFO₃ are the temperature, soil water, nitrogen, transplanting shock, phosphorous, potassium, and O₃ stress factors on photosynthesis, respectively; and SLPF is the soil fertility factor (Jones and Kiniry, 1986; Ritchie et al., 1987; Jones et al., 2003).

Leaf senescence acceleration due to O₃ stress (SLFO₃) was added to the models using Eq. (5):

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{\mathrm{SLFO}}_{\mathrm{3}}& =max\left(\mathrm{0.0},-\left({\displaystyle \frac{{\text{SFOZ}}_{\mathrm{1}}}{\mathrm{1000}}}\right)\cdot {\text{OZON}}_{\mathrm{7}}\right.\\ & \left.+\left(\mathrm{1.0}+\left({\displaystyle \frac{{\text{SFOZ}}_{\mathrm{1}}}{\mathrm{1000}}}\right)\cdot \mathrm{25.0}\right)\right),\end{array}\end{array}$$

where SFOZ₁ is the O₃ stress parameter for leaf senescence calibrated for different O₃ sensitivities of cultivars divided by a decimal correction factor of 1000 (to ensure the SFOZ₁ parameter value ranges between 0.0 and 1.0 in the model ecotype file). The SLFO₃ factor was then included in the existing daily rate of leaf area senescence function (PLAS, cm² d⁻¹) within the models as shown in Eq. (6) for maize and Eq. (7) for rice:

$$\begin{array}{}\text{(6)}& {\displaystyle }\begin{array}{rl}{\mathrm{PLAS}}_{\mathrm{maize}}& =\left(\text{PLA}-\text{SENLA}\right)\cdot \left(\mathrm{1}\right.\\ & -\min \left(\text{SLFW},\text{SLFC},\text{SLFT},\text{SLFN},\right.\\ & \left.\left.\text{SLFP},{\text{SLFO}}_{\mathrm{3}}\right)\right),\end{array}\text{(7)}& {\displaystyle }\begin{array}{rl}{\mathrm{PLAS}}_{\mathrm{rice}}& =\left(\text{PLA}-\text{SENLA}\right)\cdot \left(\mathrm{1}\right.\\ & -\min \left(\text{SLFW},\text{SLFC},\text{SLFT},\text{SLFN},\right.\\ & \left.\left.\text{SLFP},\text{SLFK},{\text{SLFO}}_{\mathrm{3}}\right)\right),\end{array}\end{array}$$

where PLA is daily plant leaf area (cm² plant⁻¹); SENLA is daily normal leaf senescence (cm² plant⁻¹); and SLFW, SLFC, SLFT, SLFN, SLFP, SLFK, and SLFO₃ are the leaf senescence stress factors due to water, light competition, temperature, nitrogen, phosphorous, potassium (CERES-Rice only), and O₃ stress, respectively (Jones and Kiniry, 1986; Ritchie et al., 1987; Jones et al., 2003).

2.2.2 CROPGRO-Soybean model

The effects of O₃ were incorporated into the CROPGRO-Soybean model using a similar approach to that described in the CERES crop models. O₃ was added into the model using the same FO₃ and PRFO₃ factors as in Eqs. (1) and (2) (for Eq. 2, PCO₂ is called PRATIO in CROPGRO-Soybean). However, CROPGRO-Soybean calculates daily photosynthesis differently than the other models do and has two different photosynthesis calculation options: leaf or canopy photosynthesis (Wilkerson et al., 1983; Boote and Pickering, 1994; Jones et al., 2003). This study focuses on the default leaf photosynthesis calculation option (which was modified to read in the CO₂ ratio effect for the PRFO₃ interaction). The daily gross photosynthesis (PG, g [CH₂O] m⁻² d⁻¹) within the model was calculated based on the limiting photosynthesis stress factors using Eq. (8) for leaf photosynthesis and Eq. (9) for canopy photosynthesis:

$$\begin{array}{}\text{(8)}& {\displaystyle }\begin{array}{rl}{\mathrm{PG}}_{\mathrm{leaf}}& =\left({\displaystyle \frac{\text{PGDAY}}{\mathrm{44.0}}}\cdot \mathrm{30.0}\cdot \text{SLPF}\right)\\ & \cdot min\left(\text{SWFAC},{\text{PRFO}}_{\mathrm{3}}\right)\cdot {\text{PSTRES}}_{\mathrm{1}},\end{array}\text{(9)}& {\displaystyle }\begin{array}{rl}{\mathrm{PG}}_{\mathrm{canopy}}& =\text{PTSMAX}\cdot \text{SLPF}\cdot {\mathrm{PG}}_{\mathrm{FAC}}\cdot {\mathrm{TPG}}_{\mathrm{FAC}}\cdot E_{\mathrm{FAC}}\\ & \cdot \text{PGSLW}\cdot \text{PRATIO}\cdot \text{PGLFMX}\\ & \cdot min\left(\text{SWFAC},{\text{PRFO}}_{\mathrm{3}}\right),\end{array}\end{array}$$

where PGDAY is daily potential photosynthesis (g [CH₂O] m⁻² d⁻¹) and SWFAC, PSTRES₁, and PRFO₃ are the soil water, phosphorous, and O₃ stress factors on photosynthesis, respectively. PTSMAX is the potential amount of CH₂O that can be produced for the full canopy (g [CH₂O] m⁻² d⁻¹), PG_FAC is a factor to compute daily PG as a function of leaf area index, TPG_FAC is a reduction factor due to less-than-optimal daytime temperature, E_FAC is the effect of nitrogen and phosphorous stress on daily canopy photosynthesis, PGSLW is the relative effect of leaf thickness on daily canopy photosynthesis, and PRATIO is the relative effect of atmospheric CO₂ on daily canopy photosynthesis (Boote and Pickering, 1994).

Leaf senescence acceleration due to O₃ stress (SLFO₃) was added to CROPGRO-Soybean using Eq. (10):

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{\text{SLFO}}_{\mathrm{3}}& =max\left(\mathrm{0.0},\left({\displaystyle \frac{{\text{SFOZ}}_{\mathrm{1}}}{\mathrm{1000}}}\right)\cdot {\text{OZON}}_{\mathrm{7}}\right.\\ & \left.-\left(\left({\displaystyle \frac{{\text{SFOZ}}_{\mathrm{1}}}{\mathrm{1000}}}\right)\cdot \mathrm{25.0}\right)\cdot \text{WTLF}\right),\end{array}\end{array}$$

where WTLF is the dry mass of leaf tissue (g_leaf m⁻²). The CROPGRO leaf senescence routine is based on existing WTLF using a different approach than the CERES leaf senescence reduction factor, so SLFO₃ has the opposite trend when compared to the CERES model calculation (Fig. 1). The SLFO₃ factor was then included in the existing daily defoliation due to daily leaf senescence (SLDOT, g m⁻² d⁻¹) calculation within the model as shown in Eq. (11):

$$\begin{array}{}\text{(11)}& \text{SLDOT}={\text{SLDOT}}_{\mathrm{n}}+max\left(\text{SLNDOT},{\text{SLFO}}_{\mathrm{3}}\right),\end{array}$$

where SLDOT_n is the natural daily leaf senescence and SLNDOT and SLFO₃ are the daily leaf senescence due to water and O₃ stress (g m⁻² d⁻¹), respectively.

2.2.3 DSSAT-NWheat model

The incorporation of O₃ into the NWheat crop model was described and validated in Guarin et al. (2019) and was used as the reference for the maize, rice, and soybean models. The approach used the same FO₃ and PRFO₃ equations as in Eqs. (1) and (2) (note that the NWheat equations were simplified from Guarin et al. (2019) by the decimal correction factor and single FOZ₁ parameter as in Eq. (1) for consistency among all models) and a similar SLFO₃ shown in Eq. (12):

$$\begin{array}{}\text{(12)}& {\mathrm{\text{SLFO}}}_{\mathrm{3}}=\left({\displaystyle \frac{{\text{SFOZ}}_{\mathrm{1}}}{\mathrm{10}}}\right)\cdot {\text{OZON}}_{\mathrm{7}}+\left(\mathrm{1.0}-\left({\displaystyle \frac{{\text{SFOZ}}_{\mathrm{1}}}{\mathrm{10}}}\cdot \mathrm{25.0}\right)\right).\end{array}$$

The O₃ effect for the different cultivar sensitivities is controlled by the FOZ₁ and SFOZ₁ parameters, as in the other models (the SFOZ₁ parameter is divided by 10 to ensure that the value ranges between 0.0 and 1.0 in the model ecotype file). The decimal correction factors vary between the crop models because the different models calculate stresses using different magnitudes.

The FO₃ and SLFO₃ responses calculated over increasing M7 O₃ concentrations are illustrated for each model in Fig. 1 using the parameter values for different O₃ cultivar classifications shown in Table 1. The FOZ₁ and SFOZ₁ parameter values for all models were determined from the cultivar sensitivities observed in the field experiments (Sect. 2.3) and the sensitivities derived from the O₃ exposure relationships from the literature (Sect. 2.5).

Figure 1 Functions for the O₃ photosynthesis reduction factor without interaction of water deficit stress and CO₂ fertilization effect (FO₃) (first column) and the O₃ leaf senescence acceleration stress factor (SLFO₃) (second column) under increasing mean 7 h (M7) O₃ concentrations for the (a, b) CERES-Maize, (c, d) CERES-Rice, (e, f) CROPGRO-Soybean, and (g, h) NWheat models. Each figure shows three different O₃ sensitivity cultivar classifications derived from the O₃ exposure–yield responses from the literature: tolerant (solid blue line), intermediate (gold short-dashed line), and sensitive (magenta long-dashed line). SLFO₃ for CROPGRO-Soybean (Eq. 10) is shown with leaf tissue dry mass (WTLF) of 1 g m⁻² for simplicity. Steeper slopes indicate a higher sensitivity to O₃ for both FO₃ and SLFO₃. Table 1 shows the parameters used in the equations for each classification of O₃ sensitivity (Eqs. 1, 5, 10, and 12).

Table 1 Summary of the O₃ photosynthesis stress parameters (FOZ₁) and the O₃ leaf senescence stress parameters (SFOZ₁) used in the FO₃ and SLFO₃ calculations (Eqs. 1, 5, 10, and 12) for the four DSSAT models under three different O₃ sensitivity cultivar classifications. The CERES and CROPGRO parameter values were determined from the O₃ exposure–yield responses in the literature (Figs. S2 and S3 in the Supplement). NWheat parameter values were from Guarin et al. (2019) and confirmed with the literature.

2.3 Observed O₃ exposure field experiments

In general, detailed field experiments of crop growth under elevated O₃ conditions for different crops are scarce and limit the granularity of model calibration. All field experiments examined in this study used dominant management conditions to limit other stresses besides O₃, e.g., water deficit or N stress, so the simulations assumed negligible outside stresses. For each crop, the DSSAT phenological and growth parameters were calibrated based on the observed control treatment with minimal O₃ stress to ensure that the models were functioning properly regardless of O₃ impact. Then, the O₃ response parameters, FOZ₁ and SFOZ₁, were calibrated based on the observed O₃ exposure–yield response between the elevated O₃ treatments and the control to simulate the O₃ effect.

For maize, the FACE experiment conducted at Champaign, Illinois, USA (40.03° N, 88.27° W; 230 m elevation), in 2018 was used for calibrating the CERES-Maize model (Choquette et al., 2020). The maize FACE experiment consisted of six cultivars grown under an ambient and an elevated O₃ treatment with n=4 (Table 2). Since there was only 1 year of data, the model was validated against the O₃ exposure-relative yield response functions from the literature (Sect. 2.5). The daily maximum temperature (TMAX), minimum temperature (TMIN), and precipitation (RAIN) weather data were collected from the nearby National Oceanic and Atmospheric Administration (NOAA) Willard Airport weather station, and the daily incoming solar radiation (SRAD) was collected from the National Aeronautics and Space Administration (NASA) Prediction Of Worldwide Energy Resources (POWER) database (https://power.larc.nasa.gov/, last access: 14 September 2022). The soil consisted of the Drummer silty clay loam soil series, and the soil parameters for this series were obtained from the United States Department of Agriculture (USDA) Natural Resources Conservation Service (NRCS) Web Soil Survey database (Table S1 in the Supplement) (NRCS, 2023). The cultivars were planted in two 3.5 m rows with a row spacing of 0.76 m on 13 May 2018 (Choquette et al., 2020). The hourly O₃ fumigation (from 10:00 to 18:00) began on 25 May 2018 and ended on 14 August 2018 and was used to calculate the daily M7 O₃ concentrations. The cultivar plots were harvested at maturity on 21 September 2018. N and water deficit stress were reported to be non-limiting, so the simulations used the non-limiting N setting within the model, and the simulated water stress was confirmed to be non-limiting with the provided rainfall. The DSSAT cultivar parameters were calibrated for phenology and growth under negligible stress conditions using the treatment with the ambient O₃ concentration (38 ppb) for each cultivar. After the phenology and growth cultivar parameters were calibrated, the FOZ₁ and SFOZ₁ O₃ response parameters were calibrated using the yield response from the elevated O₃ concentration treatments (Fig. S1a).

For soybean, data from the FACE experiment conducted at the same location in Champaign, Illinois, USA (40.03° N, 88.27° W; 230 m elevation), in 2009 and 2010 were used for model testing (Betzelberger et al., 2012). The 2009 data were used for model calibration, and the 2010 data were used for model validation. These data were previously used to incorporate O₃ effects on leaf photosynthesis into the JULES model (Leung et al., 2020). The SoyFACE experiment consisted of seven soybean cultivars grown under nine O₃ treatments with different target concentrations (Table 2). The hourly O₃ fumigation data (plots fumigated for 8 to 9 h daily except when leaves were wet) for each treatment were recorded in situ and were used to calculate the daily M7 O₃ concentrations (Betzelberger et al., 2012). The weather data were collected from the same sources as used in the maize experiment (NOAA and NASA POWER), and the soil consisted of either the Drummer silty clay loam or the Flanagan silt loam series which were obtained from the USDA NRCS Web Soil Survey database (Table S1). The initial soil conditions of the simulations were set at 95 % available water content and 100 kg N ha⁻¹ to minimize water and N stress. The cultivars were planted in plots eight rows wide and 5.4 m long, with a row spacing of 0.38 m, on 9 June 2009 and 27 May 2010. The O₃ fumigation started on 29 June 2009 and 6 June 2010 and ended on 27 September 2009 and 17 September 2010. The cultivar plots were harvested at maturity on 20 October 2009 and 30 September 2010. For each specified cultivar maturity group (Betzelberger et al., 2012), the corresponding default DSSAT maturity group parameters were used as reference and then calibrated for phenology and growth under negligible stress using the treatment with the ambient O₃ concentration (37 ppb). After the phenology and growth cultivar parameters were calibrated, the FOZ₁ and SFOZ₁ O₃ response parameters were calibrated using the yield response from the elevated O₃ concentration treatments (Fig. S1b). The parameters for both maize and soybean were calibrated using the using the one-factor-at-at-time method (Morris, 1991) until the best fit was found for the phenology, aboveground biomass and yield, and relative yield loss for each cultivar across all O₃ treatments.

Table 2 O₃ fumigation target concentration and average mean 7 h (M7, 09:00–15:59) O₃ concentrations for the 2018 maize FACE experiment (Choquette et al., 2020) and the 2009 and 2010 soybean SoyFACE experiments (Betzelberger et al., 2012).

For rice, there were no O₃ field experiment data readily available; thus a representative rice-producing location in the main North American rice-producing area in Stuttgart, Arkansas, USA (34.50° N, 91.55° W; 60 m elevation) (USDA NASS, 2010), was simulated with the default DSSAT North American rice cultivar. The year 2009 was selected for consistency with the soybean simulations. The weather data were collected from the NASA POWER database, and the dominant soil series for Arkansas County, Dewitt silt loam, was determined from the USDA NRCS Web Soil Survey database (Table S1) (NRCS, 2023). The initial soil conditions of the simulations were set at 100 % available water content and 100 kg N ha⁻¹ to ensure negligible water and N stress. A total of four 50 kg N ha⁻¹ fertilizer applications were applied throughout the season to ensure negligible N stress in the simulations. The cultivar was planted on 20 April 2009 based on the most active planting dates recorded for Arkansas in the USDA Field Crops handbook (USDA NASS, 2010), and the harvest date was automatically calculated based on when the model simulations reached physiological maturity. The default DSSAT North American rice cultivar parameters were used, and the FOZ₁ and SFOZ₁ O₃ response parameters were calibrated using the yield response from the elevated O₃ exposure functions from the literature (Sect. 2.5).

For wheat, the NWheat model was calibrated and validated using an air exclusion system O₃ exposure wheat field experiment conducted in Wake County, North Carolina, USA (35.73° N, 78.68° W; 116 m elevation), and is described in detail in Guarin et al. (2019).

2.4 Sensitivity analysis of O₃ equations and parameters

A sensitivity analysis for maize, rice, and soybean was conducted using simulations of nine constant daily M7 O₃ concentrations of 25, 40, 50, 60, 70, 80, 90, 100, and 120 ppb with different FOZ₁ and SFOZ₁ parameter values under combinations between normal or 50 % reduced rainfall and 350 or 550 ppm CO₂ concentrations to confirm that the O₃ modifications and stress interactions within the models were behaving as expected. The simulated locations and management setup for each crop were the same as the field experiments described above (Sect. 2.3). For wheat, the sensitivity analysis was based on the 1993 FACE experiment conducted in Maricopa, Arizona (33.06° N, 111.98° W; 361 m elevation) (Hunsaker et al., 1996; Kimball et al., 1999, 2017). The simulation setup for the Maricopa FACE experiment used the same 9 M7 O₃ concentrations with either a “Wet” irrigation schedule (total of 629 mm sub-surface drip irrigation at 0.23 m from planting to harvest) or a “Dry” irrigation schedule (total of 347 mm sub-surface drip irrigation at 0.23 m from planting to harvest) under 350 and 550 ppm CO₂ concentrations to examine the O₃–CO₂–water interactions as detailed in Guarin et al. (2019). For all crops, each O₃ parameter was first tested independently to examine the individual effects on photosynthesis and leaf senescence; i.e., when examining FOZ₁, SFOZ₁ was set to zero and vice versa.

2.5 Observed O₃ exposure relationships based on the literature

To confirm that the models were able to reproduce the observed relative yield loss due to O₃ stress, the simulated results were compared to well-known literature reports of O₃ exposure metrics and yield response for each crop using the M7 O₃ concentrations. The simulated locations and management conditions were the same experimental conditions as described above for each crop. For each crop, different O₃ classifications of cultivar sensitivities were defined based on more severe response to O₃ stress, i.e., tolerant, intermediate, and sensitive. These classifications of cultivar O₃ sensitivity were determined using the extensive literature review data from Mills et al. (2018b) combined with the maize and soybean FACE data for a total of 9 maize cultivars, 50 rice cultivars, 49 soybean cultivars, and 23 wheat cultivars. The literature review consisted of O₃ exposure experiments conducted in open-top chambers, experimental fields, or greenhouses and included the experiments that contributed to the widely applied Weibull O₃ response function (Heck et al., 1984; Adams et al., 1989; Lesser et al., 1990; Wang and Mauzerall, 2004; Tai et al., 2021; Feng et al., 2022). The selection criteria of the data are described in detail in Mills et al. (2018b).

The yield data from the literature experiments were standardized as performed by Mills et al. (2018b) and described by Osborne et al. (2016). For each experiment, linear regression was used to determine the yield at 25 ppb M7 O₃, and this value was the reference for calculating the relative yield; i.e., relative yield was calculated as the actual observed yield divided by the yield at 25 ppb O₃. The 25 ppb M7 O₃ threshold was chosen for proper comparison to the model results. After calculating the yield relative to 25 ppb M7 O₃, a linear regression for each cultivar was performed using R statistical software, v4.3.0 (R Core Team, 2023; Wickham, 2016; Wickham et al., 2023), to determine the O₃ exposure response (Fig. S2). The cultivar O₃ exposure responses were then classified into three evenly distributed quantiles, 0 %–33 %, 33 %–66 %, and 66 %–100 %, chosen to represent the three O₃ sensitivity classifications, sensitive, intermediate, and tolerant, respectively (Fig. S3). These data were used to determine the model FOZ₁ and SFOZ₁ values of each of the O₃ cultivar classifications shown in Table 1 to evaluate if the models could accurately reproduce the O₃ exposure–yield responses.

3 Results

3.1 Calibration of crop models and simulated relative yield loss against O₃ exposure field experiments

The simulated phenology (anthesis (flowering) and physiological maturity dates), biomass, yield, and relative yield due to elevated O₃ stress from the maize and soybean experiments were compared to the field observations to determine performance of the O₃ equations within the models (Tables 3–5; Figs. 2, 3, and S1). The relative yield due to O₃ stress was calculated by dividing the yield of each corresponding O₃ treatment over the control yield, i.e., the baseline O₃ treatment, and multiplying by 100 to convert to a percentage. The relative yield loss was the difference between 100 % and the calculated relative yield. There were no O₃ field experiment data for rice, so the rice O₃ parameter values and performance were compared to the O₃ exposure–yield response functions from the literature (Sect. 3.3).

The maize and soybean cultivars had different sensitivities to O₃ stress which were accounted for by using different FOZ₁ and SFOZ₁ values (Fig. S1). The calibrated CERES-Maize and CROPGRO-Soybean models simulated the physiological maturity within 4 d of the observations (Table 5; root-mean-square error (RMSE) = 0.0 d for maize 2018, 3.70 d for soybean 2009, and 3.30 d for soybean 2010). The calibrated CERES-Maize model was able to reproduce the yield and relative yield loss very well across all six cultivars (Fig. 2; RMSE = 107 kg ha⁻¹ and 2 %; r²=0.99 and 0.99, respectively). This ideal model performance was because only two O₃ treatments were available for each maize cultivar which simplified the calibration process (Fig. S1a). The CROPGRO-Soybean model was able to reproduce the biomass, yield, and relative yield loss due to O₃ stress well for the calibration year, 2009 (Fig. 3a, b, c; RMSE = 1179 kg ha⁻¹, 328 kg ha⁻¹, and 10 %; r²=0.81, 0.88, and 0.85), and acceptably for the evaluation year, 2010, across all seven cultivars (Fig. 3d, e, f; RMSE = 3339 kg ha⁻¹, 1291 kg ha⁻¹, and 16 %; r²=0.59, 0.71, and 0.66). The model overestimated biomass and yield for all cultivars and treatments in 2010, which was likely the result of a factor outside of the model setup that mitigated the increased incoming solar radiation when compared to 2009 (Sect. 4.3). The calibration and evaluation for the NWheat model was conducted and validated in Guarin et al. (2019), where the model reproduced the observed relative yield due to O₃ stress with a normalized root-mean-square error (NRMSE) of 23 % and an r² of 0.94, 0.91, and 0.88 for the tolerant, intermediate, and sensitive O₃ sensitive cultivar classifications.

Table 3 CERES-Maize cultivar and O₃ parameters used to simulate the six maize cultivars from the 2018 FACE field experiment (Choquette et al., 2020). P1: thermal time from seedling emergence to the end of the juvenile phase (expressed in degree days above a base temperature of 8 °C). P2: extent to which daily development is delayed for each hour increase in photoperiod above the longest photoperiod at which development proceeds at a maximum rate (which is considered to be 12.5 h). P5: thermal time from silking to physiological maturity (expressed in degree days above a base temperature of 8 °C). G2: maximum possible number of kernels per plant. G3: kernel filling rate during the linear grain filling stage and under optimum conditions (mg d⁻¹). PHINT: phyllochron interval, i.e., the interval in thermal time (degree days) between successive leaf tip appearances. FOZ₁: O₃ effect on photosynthesis. SFOZ₁: O₃ effect on leaf senescence.

Table 4 CROPGRO-Soybean cultivar and O₃ parameters used to simulate the seven soybean cultivars based on the maturity groups defined in the SoyFACE field experiment (Betzelberger et al., 2012). CSDL: critical short-day length below which reproductive development progresses with no day length effect (for short-day plants) (hour). PPSEN: slope of the relative response of development to photoperiod with time (positive for short-day plants) (per hour). EM-FL: time between plant emergence and flower appearance (R1) (photothermal days). FL-SH: time between first flower and first pod (R3) (photothermal days). FL-SD: time between first flower and first seed (R5) (photothermal days). SD-PM: time between first seed (R5) and physiological maturity (R7) (photothermal days). FL-LF: time between first flower (R1) and end of leaf expansion (photothermal days). LFMAX: maximum leaf photosynthesis rate at 30 °C, 350 vpm CO₂, and high light (mg CO₂ m⁻² s⁻¹). SLAVR: specific leaf area of cultivar under standard growth conditions (cm² g⁻¹). SIZLF: maximum size of full leaf (three leaflets) (cm²). XFRT: maximum fraction of daily growth that is partitioned to seed and shell. WTPSD: maximum weight per seed (g). SFDUR: seed filling duration for pod cohort at standard growth conditions (photothermal days). SDPDV: average seed per pod under standard growing conditions (number per pod). PODUR: time required for cultivar to reach final pod load under optimal conditions (photothermal days). THRSH: threshing percentage. The maximum ratio of (seed per (seed + shell)) at maturity. SDPRO: fraction protein in seeds (g(protein) per g(seed)). SDLIP: fraction oil in seeds (g(oil) per g(seed)). FOZ₁: O₃ effect on photosynthesis. SFOZ₁: O₃ effect on leaf senescence.

Table 5 Observed and simulated anthesis day and maturity day for the six maize cultivars from the 2018 FACE experiment (Choquette et al., 2020) and the seven soybean cultivars from the 2009 and 2010 soybean SoyFACE experiments (Betzelberger et al., 2012). The observed maturity dates were estimated from the single reported harvest date for all cultivars, but there may have been minor variation between the different cultivars. Observed anthesis was not available for soybean.

Figure 2 CERES-Maize model calibration of the 2018 FACE O₃ field experiment conducted in Champaign, Illinois, USA (Choquette et al., 2020). Simulated and observed (a) yield and (b) relative yield due to elevated O₃ stress (compared to the ambient control treatment) for six maize cultivars (colored points). The root-mean-square error (RMSE) and coefficient of determination (r²) show the model performance across all cultivars. The solid black line shows 1:1 comparison, and the dotted black line shows linear fit across all cultivars. For maize, only 1 year of experimental data were available for calibration and evaluation. The model cultivar parameters are shown in Table 3.

Figure 3 CROPGRO-Soybean model performance and evaluation of the SoyFACE O₃ field experiment conducted in Champaign, Illinois, USA (Betzelberger et al., 2012). Simulated and observed (a, d) above-ground biomass, (b, e) yield, and (c, f) relative yield in response to the nine progressive O₃-increasing treatments (Table 2) for seven soybean cultivars (colored points). Relative yield is compared to the ambient control treatment within each year. The 2009 SoyFACE field experiment was used for model calibration (a, b, c), and the 2010 SoyFACE field experiment was used for model evaluation (d, e, f). The root-mean-square error (RMSE) and coefficient of determination (r²) show the model performance across all cultivars. The solid black line shows 1:1 comparison, and the dotted black line shows linear fit across all cultivars. The model cultivar parameters are shown in Table 4.

3.2 Sensitivity analysis and combined effects of O₃, CO₂, and water deficit stress on yields

The simulated relative yield losses due to O₃ stress increased for all crops as the M7 O₃ concentrations increased above the 25 ppb threshold when examining the photosynthesis and leaf senescence responses independently, as expected (Figs. 4–7). The simulated actual yields for all crops are shown in Tables S2–S9. Wheat was the most sensitive crop to O₃ stress of the four crops examined (compare slopes in Figs. 4–7a and b), which agrees with previous literature (Mills et al., 2018b). For each model, simulations using an FOZ₁ or SFOZ₁ example value of 0.5 were examined in more detail to illustrate the O₃–CO₂–water interactions (Figs. 4–7c and d, respectively). For all crops, the Dry/reduced rainfall and low CO₂ treatment produced the lowest yields, while the Wet/normal rainfall and high CO₂ produced the highest yields (Tables S2–S9). The simulated O₃ effect was larger when water was non-limiting; i.e., the higher rainfall and irrigated treatments experienced larger losses due to O₃ stress because of increased stomatal uptake. The simulated O₃ effect was reduced under the higher CO₂ concentrations, thus capturing the responses from stomatal closure and the photosynthetic benefits from the CO₂ fertilization effect.

Figure 4 Sensitivity analysis using the CERES-Maize model to simulate relative yield due to elevated O₃ stress (relative to 25 ppb M7 O₃) for a range of (a) the photosynthesis O₃ stress parameter (FOZ₁) and (b) the leaf senescence O₃ stress parameter (SFOZ₁) values under the normal rainfall and 350 ppm CO₂ scenario and an example of (c) FOZ₁ and (d) SFOZ₁ set at 0.5 under the 50 % reduced rainfall and 350 ppm CO₂ (solid red line), normal rainfall and 350 ppm CO₂ (solid black line), 50 % less rainfall and 550 ppm CO₂ (dashed red line), and normal rainfall and 550 CO₂ (dashed black line) scenarios. The Champaign, Illinois, USA, 2018 FACE weather, soil, and dominant management conditions were used for the reference location. Each O₃ parameter was tested independently; i.e., when examining FOZ₁, SFOZ₁ was set to zero and vice versa. The simulated actual yields are shown in Tables S2 and S3.

Figure 5 Sensitivity analysis using the CERES-Rice model to simulate relative yield due to elevated O₃ stress for a range of (a) FOZ₁ and (b) SFOZ₁ values under the normal rainfall and 350 ppm CO₂ scenario and an example of (c) FOZ₁ and (d) SFOZ₁ set at 0.5 under the 50 % reduced rainfall and 350 ppm CO₂ (solid red line), normal rainfall and 350 ppm CO₂ (solid black line), 50 % less rainfall and 550 ppm CO₂ (dashed red line), and normal rainfall and 550 CO₂ (dashed black line) scenarios. The Stuttgart, Arkansas, USA, 2009 weather, soil, and dominant management conditions were used for the reference location. Each O₃ parameter was tested independently; i.e., when examining FOZ₁, SFOZ₁ was set to zero and vice versa. The simulated actual yields are shown in Tables S4 and S5.

Figure 6 Sensitivity analysis using the CROPGRO-Soybean model to simulate relative yield due to elevated O₃ stress for a range of (a) FOZ₁ and (b) SFOZ₁ values under the normal rainfall and 350 ppm CO₂ scenario and an example of (c) FOZ₁ and (d) SFOZ₁ set at 0.5 under the 50 % reduced rainfall and 350 ppm CO₂ (solid red line), normal rainfall and 350 ppm CO₂ (solid black line), 50 % less rainfall and 550 ppm CO₂ (dashed red line), and normal rainfall and 550 CO₂ (dashed black line) scenarios. The Champaign, Illinois, USA, 2009 SoyFACE weather, soil, and dominant management conditions were used for the reference location. Each O₃ parameter was tested independently; i.e., when examining FOZ₁, SFOZ₁ was set to zero and vice versa. The simulated actual yields are shown in Tables S6 and S7. Figure S6 shows the relative biomass loss corresponding to SFOZ₁ (d) to explain the inverted CO₂ effect under the 50 % rainfall treatment.

Figure 7 Sensitivity analysis using the NWheat model to simulate relative yield due to elevated O₃ stress for a range of (a) FOZ₁ and (b) SFOZ₁ values under the Wet irrigation and 350 ppm CO₂ scenario and an example of (c) FOZ₁ and (d) SFOZ₁ set at 0.5 under the Dry irrigation and 350 ppm CO₂ (solid red line), Wet irrigation and 350 ppm CO₂ (solid black line), Dry irrigation and 550 ppm CO₂ (dashed red line), and Wet irrigation and 550 CO₂ (dashed black line) scenarios. The Maricopa, Arizona, USA, 1993 FACE weather, soil, and management conditions were used for the reference location (Kimball et al., 1999; Guarin et al., 2019). Each O₃ parameter was tested independently; i.e., when examining FOZ₁, SFOZ₁ was set to zero and vice versa. The simulated actual yields are shown in Tables S8 and S9.

3.3 Simulated relative yield loss compared to O₃ relationships in the literature

For all crops, the literature showed a large range of relative yield losses due to O₃ stress caused by different cultivar O₃ sensitivities (Fig. S2). Wheat was the most sensitive crop to O₃ stress with an average yield loss of 0.70±0.39 % (mean ± SD) per ppb M7 O₃ increase above 25 ppb, followed by soybean, maize, and then rice (average yield losses of 0.60±0.39 %, 0.39±0.26 %, and 0.32±0.37 % per ppb M7 O₃ increase above 25 ppb, respectively) (average of slopes in Table S10). To encompass the high variability of yield losses, the cultivars were classified into the O₃-tolerant, -intermediate, and -sensitive cultivar O₃ sensitivities (Fig. S3). Since the cultivar sensitivities were not originally specified in the literature, the FOZ₁ and SFOZ₁ parameters used in the models were adjusted to provide the best fit across the O₃ exposure responses (Table 1). Overall, the models reproduced the simulated O₃ exposure relationships from the literature well: the RMSE for maize, rice, soybean, and wheat across all three O₃ exposure sensitivities were 6.6 %, 7.8 %, 4.0 %, and 5.4 %, respectively (Fig. 8). The models performed better (lower RMSE) for the O₃-tolerant and O₃-intermediate cultivar sensitivities compared to the O₃-sensitive cultivar sensitivity, but all models explained the variance well (r²>0.96 across all O₃ sensitivities). This suggests that different combinations of FOZ₁ and SFOZ₁ can be calibrated for specific observations to emulate the variation in different O₃ exposure responses.

Figure 8 Simulated relative yield due to O₃ stress (solid lines) compared to the O₃ exposure relationships (dotted lines) from the literature data (symbols) for the (a) CERES-Maize, (b) CERES-Rice, (c) CROPGRO-Soybean, and (d) NWheat models. The calculated relative yield from the well-known Weibull O₃ response functions (dashed black lines, equations listed in Table S11) is based on the US NCLAN network O₃ exposure field experiments conducted between the 1960s and the 1980s (Adams et al., 1989; Lesser et al., 1990; Wang and Mauzerall, 2004; Tai et al., 2021). The O₃ exposure–yield response linear fits of the three O₃ sensitivities, tolerant (blue), intermediate (gold), and sensitive (magenta), are given in Fig. S3. The cultivars were classified by grouping the cultivar O₃ exposure–yield response (Fig. S2) into three evenly distributed quantiles: 66 %–100 %, 33 %–66 %, and 0 %–33 %, respectively. The O₃ sensitivities determined for each cultivar are listed in Table S10. The simulated results for the crop models use the FOZ₁ and SFOZ₁ values from Table 1. For each model, the same weather, soil, and dominant management conditions as in the normal rainfall and 350 ppm CO₂ treatment of the sensitivity analysis were used as reference (the O₃ response functions from the literature included O₃ field experiments conducted when the atmospheric CO₂ concentration was ∼350 ppm). The literature data consist of the relative yields (scaled to 25 ppb M7 O₃) of the cultivars examined in the Mills et al. (2018b) literature review combined with the maize and soybean cultivars used in this study for a total of 9 maize cultivars, 50 rice cultivars, 49 soybean cultivars, and 23 wheat cultivars (listed in Table S10). For the O₃ sensitivity of each crop, the root-mean-square error (RMSE) and coefficient of determination (r²) show the crop model performance (RMSE_sim) and the Weibull response function performance (RMSE_weibull) compared to the linear fit of the O₃ exposure literature data (text color corresponds to O₃ sensitivity). The color-shaded area shows the standard error for the linear fit of the literature data for each of the cultivar O₃ sensitivities.

4 Discussion

4.1 Simulating O₃ damage on crop yields

The measured yield losses for the maize FACE experiment were between 5 % and 40 % for the M7 O₃ concentrations when increasing from the ambient concentration (38 ppb) to the elevated O₃ treatment (77 ppb), a yield loss of 0.14 % to 1.01 % per ppb M7 O₃ above the ambient concentration, depending on the O₃ cultivar sensitivity (Fig. 2b). `NC338xHp301' and `Mo17xHp301' were classified as O₃-tolerant because of relatively small yield losses of 5 % and 6 %, respectively; `B73xMo17' was classified as O<sub>3</sub>-intermediate with a yield loss of 11 %; and `B73xHp301', `Mo17xNC338', and `B73xNC338' were sensitive to O₃ effects with yield losses of 22 %, 30 %, and 40 %, respectively (Fig. S1a, Table S10). These cultivar O₃ sensitivities are based on a single experimental year, so additional testing is needed to further corroborate the classifications. Overall, the calibrated CERES-Maize model was able to reproduce these observed yield losses within 1 %, i.e., simulated yield losses between 5 % and 41 %, or 0.12 % and 1.05 % per ppb O₃ increase above the ambient concentration. These yield losses were also calculated relative to 25 ppb (as described in Sect. 2.5) for consistency with the literature, which resulted in simulated yield losses between 0.12 % and 0.93 % per ppb M7 O₃ increase above 25 ppb across the six cultivars.

When comparing the simulations to the maize O₃ exposure–yield relationships from the literature, the model simulated average yield losses of 0.16 %, 0.36 %, and 0.82 % per ppb M7 O₃ increase above 25 ppb for the O₃-tolerant, -intermediate, and -sensitive cultivar O₃ sensitivities, respectively (Fig. 8a solid lines). This agreed well with the literature yield losses of 0.24 %, 0.33 %, and 0.71 % per ppb M7 O₃ increase above 25 ppb for the O₃-tolerant, -intermediate, and -sensitive cultivar sensitivities, respectively (Figs. S3a, 8a dotted lines). The O₃ parameter values used for the literature comparison were determined to provide the best fit across the literature experiments consisting of nine maize cultivars, but these O₃ parameter values could be calibrated for other scenarios and cases, i.e., higher or lower cultivar O₃ sensitivity.

The measured yield losses for the SoyFACE experiment were between 51 % and 77 % for the M7 O₃ concentrations when increasing from the ambient concentration (37 ppb) to the highest O₃ treatment (126 ppb) in 2009, a yield loss of 0.57 % to 0.86 % per ppb M7 O₃ above the ambient concentration, depending on the cultivar O₃ sensitivity (Fig. 3c). The calibrated CROPGRO-Soybean model reproduced observed yield losses within 10 %, i.e., simulated yield losses between 59 % to 80 %, or 0.66 % to 0.90 % per ppb O₃ increase. Based on the calculated O₃ classifications from the literature and low yield divergence across the seven cultivars (Fig. S1b), `Pioneer93B15', `Dwight', ` IA-3010', and `LN97-15076' were considered O₃-intermediate, and `HS93-4118', `Loda', and `Pana' were considered O₃-sensitive (Table S10). In 2010, the observed soybean yield losses ranged between 31 % and 76 % when increasing from the ambient concentration (37 ppb) to the highest O₃ treatment (84 ppb), a yield loss of 0.65 % to 1.60 % per ppb M7 O₃ above the ambient concentration. The model underestimated yield losses in 2010, between 27 % and 44 %, but, because the experimental setup was the same for both years, an external factor that was not considered in the simulations may have affected yields (Sect. 4.3). The 2010 yield losses were of a similar magnitude to the 2009 yield losses, but the 2010 experiment had higher yield loss and variation per ppb O₃ increase with lower average M7 O₃ concentrations (Table 2, Fig. S4a).

When comparing the simulations to the soybean O₃ exposure–yield relationships from the literature (Fig. 8c), an average yield loss of 0.36 %, 0.64 %, and 0.96 % per ppb M7 O₃ increase above 25 ppb was simulated for the O₃-tolerant, -intermediate, and -sensitive cultivar O₃ sensitivities, respectively. This was substantiated by the literature yield losses of 0.45 %, 0.63 %, and 0.84 % per ppb M7 O₃ increase above 25 ppb for the O₃-tolerant, -intermediate, and -sensitive cultivar O₃ sensitivities, respectively (Figs. S3c, 8c dotted lines). The literature data consisted of 49 soybean cultivars, which had a smaller range of O₃ sensitivities compared to the other crops, although there were outliers where yield increased under higher O₃ concentrations (described in Sect. 4.2).

The CERES-Rice model simulated an average yield loss of 0.05 %, 0.23 %, and 0.66 % per ppb M7 O₃ increase above 25 ppb for the O₃-tolerant, -intermediate, and -sensitive cultivar O₃ sensitivities, respectively (Fig. 8b solid lines). The rice literature had the most cultivars (50) of the four crops examined, and the simulated yield losses for the O₃-tolerant and -intermediate cultivar O₃ sensitivities agreed well with the literature yield losses of 0.07 % and 0.24 % per ppb M7 O₃ increase above 25 ppb, respectively (Fig. 8b dotted lines). A larger discrepancy between the simulated yield loss for the O₃-sensitive classification and the literature O₃-sensitive yield loss of 0.49 % per ppb M7 O₃ increase above 25 ppb was due to the higher variability within the literature data (Fig. 8b shaded area).

Using the calibrated NWheat model, the simulated yield losses were 0.26 %, 0.66 %, and 1.23 % per ppb M7 O₃ increase above 25 ppb for the O₃-tolerant, -intermediate, and -sensitive cultivar O₃ sensitivities, respectively (Fig. 8d). These simulated yield losses were corroborated by the reported average yield losses of 0.33 %, 0.61 %, and 1.11 % per ppb M7 O₃ increase above 25 ppb for the O₃-tolerant, -intermediate, and -sensitive cultivar O₃ sensitivities, respectively. The literature expanded across different ranges of O₃ concentrations for all crops, and yield loss per ppb is not always constant over an expansive range of O₃ concentrations, so the model O₃ parameter values can be adjusted for higher or lower cultivar O₃ sensitivity.

As an additional check of model performance, the calculated relative yield from the well-known Weibull O₃ response functions (Table S11) were compared to the literature O₃ exposure linear yield responses for each crop and O₃ classification (Fig. 8). The Weibull function performance was then compared to the simulated crop model results. Overall, the crop model simulations performed better (lower RMSE and higher r²) than the Weibull response functions across all crops for all three O₃ classifications, except the O₃-intermediate classification for soybean which had <1 % difference between the RMSE (compare RMSE and r² in Fig. 8). The performance results suggest that it is best to use calibrated crop models when available and that the Weibull response functions are mainly representative of O₃-intermediate classifications for maize, rice, and soybean and O₃-tolerant classifications for wheat.

4.2 Simulated relative yield loss with the combined effects of O₃, CO₂, and water deficit stress

The sensitivity analyses showed that the yield losses due to O₃ stress were higher under the normal rainfall and low CO₂ treatment, which agrees with previous literature that increased water availability increases O₃ impact due to increased stomatal uptake (Khan and Soja, 2003; Biswas et al., 2013). It was unexpected that the simulated O₃ photosynthetic response difference between the normal and reduced rainfall treatments for maize was less than 1 % (Fig. 4c). This was because the model simulated low water deficit stress under the 50 % reduced rainfall treatment which obscured the O₃–water stress dynamics. Further reducing the rainfall to 40 % of the normal amount increased the simulated water deficit stress and produced the photosynthetic O₃–water dynamics consistent with the other models (Fig. S5). The elevated CO₂ concentration mitigated the detrimental effect of O₃ stress in the photosynthetic response for all models (Figs. 4–7c), which agrees with recent global findings that elevated CO₂ concentrations can mitigate and even negate elevated O₃ impacts (Xia et al., 2021; Tai et al., 2021). Interestingly, the CROPGRO-Soybean model simulated an inverse O₃–CO₂ effect on relative yield under the 50 % rainfall condition when examining SFOZ₁ in detail (Fig. 6d). This inverse yield response was due to the low actual yield simulated under the 50 % rainfall and low CO₂ treatment (<2000 kg ha⁻¹, Table S7) which resulted in smaller changes in yield compared to the 50 % rainfall and high CO₂ treatment, but the overall simulated aboveground biomass O₃–CO₂–water interaction was as expected (Fig. S6).

For several of the observations from the actual soybean field experiment using `Pana', the yield increased under higher O<sub>3</sub> concentrations (∼2 % to 18 %; Figs. 3c, S1b, and S7). In some cases it is possible that elevated O<sub>3</sub> concentrations can benefit a crop via hormesis, a process where low levels of intermittent stress may benefit overall crop growth through improved resiliency (Calabrese, 2014). It is also possible that if elevated O<sub>3</sub> concentrations reduce biomass growth throughout the season, and therefore reduce nutrient resource demand throughout the season, small yield increases can occur from a larger pool of resources available during the key reproductive/grain filling period (Asseng and Van Herwaarden, 2003; Guarin et al., 2019). This increase in yield under higher O<sub>3</sub> concentrations was also observed under several other soybean and rice cultivars from the literature (Fig. S2b and c). However, a soybean cultivar from the literature, `Cumberland', was reported to have a 34 % increase under elevated O₃ (67 ppb) compared to the control treatment (25 ppb), but such a large increase may indicate that another outside factor affected the yields. Mulchi et al. (1988) speculated that the large yield difference was due to changes in the seasonal water dynamics, thereby causing increased drought stress under the control treatment compared to the elevated O₃ treatment. Reproducing rare occurrences where elevated O₃ may result in yield increases can be a challenge for the models because of the linear response of the stress equations (Fig. S7). However, it may be possible depending on the simulated interactions between seasonal dynamics of resources, as shown with the sensitivity analysis of wheat yields in Guarin et al. (2019).

4.3 Uncertainty in model simulations and O₃ exposure field experiments

Crop models contain uncertainties due to simplification of complex biological processes, but field experiments may also contribute uncertainty via measurement. The soybean simulations overestimated both biomass and yield across all cultivars and treatments for the 2010 SoyFACE experiment. Since both the ambient and elevated O₃ treatments were overestimated, it is unlikely that the simulated O₃ interactions caused the discrepancy. Examining the weather input showed a 14 % increase in cumulative incoming solar radiation for the 2010 growing season compared to the 2009 growing season (Fig. S4b). The 2010 season was warmer than the 2009 season, average seasonal temperature of 23.4 °C compared to 19.1 °C, but no heat stress was reported, and the difference in rainfall was negligible, 445 mm compared to 454 mm. Since management was the same for both years, and no water or N stresses were reported, it was expected that the 2010 yields would be higher than the 2009 yields due to the increased solar radiation, but the average 2010 yield across all cultivars for the ambient treatment decreased, 3300 kg ha⁻¹ in 2010 compared to 3700 kg ha⁻¹ in 2009. Therefore, it is possible that an outside stress factor not considered within the model limited soybean growth in the field in 2010 which led to the model overestimating biomass and yield. One possibility is that increased rainfall during the beginning of the 2010 season (221 mm in first 30 d compared to 153 mm in first 30 d of the 2009 season; Fig. S4c) may have resulted in germination or emergence stress due to excessive water such as flooding or lodging, which are factors not yet considered in the crop models.

The sensitivity analyses showed that the CO₂ effect was more pronounced in the model photosynthesis response than in the leaf senescence response (compare solid and dashed lines in Figs. 4–7c and d). This is because the models do not have a CO₂ effect directly applied to the daily leaf senescence calculation, whereas CO₂ directly affects the daily photosynthesis calculation (PCARB in Eqs. 3 and 4 and PRATIO in Eq. 9). Improved CO₂ representation within the crop models is being explored through the Agricultural Model Intercomparison and Improvement Project (AgMIP) studies (Ahmed et al., 2017, 2019; Toreti et al., 2020), but additional high-quality data is needed for model testing.

5 Conclusion

Crop responses to elevated O₃ concentrations were incorporated into the DSSAT CERES-Maize, CERES-Rice, CROPGRO-Soybean, and NWheat crop models via functions reducing photosynthetic activity and accelerating leaf senescence. Model testing showed that each of the four models reproduced the observed O₃ response from field experiments and previous literature and the expected interactions between O₃, CO₂, and water deficit stress. The simulated yield responses were also more representative of the O₃ exposure literature data than the well-known Weibull O₃ response functions for all crops. Thus, this incorporation allows improved simulation of the heterogeneity of O₃ impacts across geographical regions and systems, and across years within seasons, which is more representative of real-world interactions than using a generic damage coefficient. Overall, increasing M7 O₃ concentrations had a negative effect on growth and yield across all four crops, and this negative effect was exacerbated by increased water availability and ameliorated by elevated CO₂ concentrations. The O₃ impact and stress response of the crop depends on the stress severity, duration, frequency, cultivar sensitivity, and seasonal timing (i.e., developmental stage), which can be accounted for by using the updated crop models.

The addition of O₃ stress functionality into crop models will improve both near- and long-term simulations of global environmental interactions using a key factor that is often not included in agricultural and climate change assessments. The DSSAT models in this study can be used to simulate the O₃ impacts on crops in combination with climate change. The O₃ parameter values in this study can be used as preliminary approximations, but, to further improve model performance and robustness of the O₃ stress routines, the models and parameters should continue to be tested and calibrated with additional O₃ exposure experimental data when available. In addition, the models should be compared with other O₃-modified crop models as part of multi-model ensemble intercomparison and improvement assessments conducted by the AgMIP (https://agmip.org/, last access: 28 March 2024). As a next step, the AgMIP ozone team is currently conducting a multi-model ensemble study with crop models that have the capacity to evaluate the responses of future crop yields to different ozone concentrations. This effort will help produce more robust estimates of climate change impacts in global agriculture. The framework described here can be used by other process-based crop models, local or gridded, to incorporate O₃ stress interactions into the model. This model improvement also suggests potential future collaboration between crop modelers and remote sensing experts using weather and climate models with dynamic chemistry components, such as the NASA Atmosphere Observing System (https://aos.gsfc.nasa.gov/, last access: 28 March 2024).