A first calibration of JULES-crop version 7.4 for rice using the novel O₃-FACE experiment in China

Ozone (O₃) pollution poses an escalating threat to rice production and food security in China, with concentrations projected to rise under future climate scenarios. Accurately quantifying O₃ impacts on rice is thus crucial for informed agricultural planning. This study is the first to utilise free-air concentration enrichment (FACE) observations specific to rice for calibrating a crop model (Joint UK Land Environment Simulator with crops, JULES-crop) and assessing the impacts of O₃. FACE experiments, which involve growing crops under natural field conditions while exposing them to elevated O₃ levels, provide an ideal approach for studying the effects of O₃ on crops. Utilising data from the only O₃-FACE facility dedicated to rice, we calibrated physiological and O₃-response parameters in JULES-crop and evaluated the model against additional independent FACE observations. The calibration establishes this as the first crop model refined with ideal open-air field observations, significantly enhancing its capability to simulate rice growth processes and O₃-induced yield losses and surpassing the performance of simulations based on the default parameters in JULES-crop. With this newly calibrated model, JULES-crop is now equipped to assess the impacts of O₃ on agriculture, offering a valuable tool to inform mitigation strategies.

1 Introduction

Rice is a staple food for over half of the world's population and plays a crucial role in global food security. The rising concentration of ozone (O₃) is a major concern, contributing to significant losses in crop production worldwide (Van Dingenen et al., 2009). Mills et al. (2018) estimated that the average global yield loss of rice due to O₃ was 4.4 % between 2010 and 2012. In China, O₃ caused relative rice yield losses of 6.2 %–52.9 % between 2014 and 2018 and of 23 % between 2017 and 2019 (Feng et al., 2022; Xu et al., 2021). Consequently, assessing the impact of O₃ on rice growth is essential, especially as O₃-polluted areas overlap with crop-growing regions and pose a long-term threat to food security (Emberson et al., 2018).

The main O₃ dose-response functions used to assess rice yield loss include concentration-based methods, such as the accumulated dose of O₃ over 40 ppb (AOT40) and the daily mean 7 h concentrations (M7), and flux-based methods, such as the phytotoxic O₃ dose (POD) (Tai et al., 2021). Both concentration-based and flux-based methods can establish a relationship with relative yield loss based on field experiments. The relationship between relative yield loss and O₃ level, known as the O₃ response function, is a valuable tool that underpins extensive research into crop yield losses caused by O₃ exposure (Ramya et al., 2023).

Some crop models have incorporated O₃ parameters to better understand O₃'s impacts (Guarin et al., 2024; Leung et al., 2020; Ewert and Porter, 2000). For instance, the Decision Support System for Agrotechnology Transfer (DSSAT) crop model established an O₃ stress factor using the M7 metric (Guarin et al., 2024). GLAM-ROC simulated O₃ effects by reducing evapotranspiration, transpiration efficiency, and the harvest index based on AOT40 metric (Droutsas et al., 2020). The Joint UK Land Environment Simulator with crops (JULES-crop) integrated a flux-based O₃ damage scheme developed by Sitch et al. (2007) to assess reductions in net photosynthesis. Flux-based methods account for stomatal conductance and environmental conditions, such as temperature and the vapour pressure deficit, to modify O₃ uptake and thus directly link the absorbed O₃ dose to physiological damage. Compared with concentration-based methods, flux-based methods exhibit enhanced performance in correlating O₃ levels with relative yield loss, enabling more precise assessments (Pleijel et al., 2004, 2022; Mills et al., 2011; Ronan et al., 2020). Nonetheless, O₃-related parameters in crop models require calibration to ensure reliable performance, even when using a flux-based O₃ scheme.

Open-top chambers (OTCs) and free-air concentration enrichment (FACE) experiments are two major methods used to help calibrate parameters in crop models. State-of-the-art FACE experiments, which provide more natural environments for crops, are ideal for establishing O₃ exposure metrics and investigating the impacts of O₃ on crops (Montes et al., 2022; Feng et al., 2018). To date, only four O₃-FACE facilities have been established for crops worldwide (Montes et al., 2022): wheat and rice experiments in China (Tang et al., 2011), wheat experiments in India (Yadav et al., 2019), grape experiments in Italy (Moura et al., 2023), and soybean experiments in the United States (Aspray et al., 2023). However, the rice-specific O₃-FACE experiment has not yet been used to calibrate any crop models.

The parameterisation of crops in JULES was developed by Osborne et al. (2015). JULES-crop incorporates flux-based O₃ exposure metrics to analyse the loss of accumulated carbon based on the exact O₃ flux entering the crop stomata, which is influenced by environmental conditions (Sitch et al., 2007). The impact of O₃ on crops is also reflected in reductions in crop height, leaf area index (LAI), and crop yields. Additionally, Tai et al. (2021) highlighted that mechanistic crop models such as JULES-crop can combine the fertilisation effects of atmospheric carbon dioxide (CO₂) with O₃ influence. Thus, JULES-crop is a suitable tool for investigating the effects of O₃ on crops, accounting for environmental factors that modify the mechanisms of O₃ effects (Leung et al., 2022). However, the crop growth and development parameters for rice, as well as the O₃ impact parameters within JULES-crop, have not yet been calibrated. Calibrating JULES-crop would enhance its performance in simulating rice production under O₃ influence.

In this research, we calibrated the rice parameters in JULES-crop using novel O₃-FACE data, enabling leading-edge future assessments of O₃ damage to rice. The study has three key objectives: (1) to calibrate JULES-crop using novel O₃-FACE field data; (2) to evaluate the model's performance in capturing crop growth characteristics using independent observations; and (3) to assess the impact of O₃ on rice physiology, phenology, and yields. This research enhances understanding of the mechanisms through which O₃ affects rice growth and development, providing a stronger basis for characterising the future impact of O₃ on rice production.

2 Method

2.1 Description of JULES-crop

JULES-crop is an extension of JULES, a land surface model designed to simulate the fluxes of carbon, water, energy, and momentum between the land surface and the atmosphere (Best et al., 2011; Clark et al., 2011). JULES-crop was developed to simulate the growth and development of major crops, including wheat, soybean, maize, and rice, under a range of environmental influences such as temperature, precipitation, radiation, and soil moisture (Osborne et al., 2015). Its structure, illustrated in Fig. 1, incorporates the physiological processes of crops, including photosynthesis, respiration, and biomass accumulation.

Figure 1 Schematic of JULES-crop.

JULES-crop simulates the physiological and phenological processes of crops, predicting yields at both field and global scales. This capability makes it a valuable tool for understanding the impacts of climate change and air pollution on agriculture (Leung et al., 2022; Wolffe et al., 2021; Vianna et al., 2022). To date, winter wheat (in preparation), maize (Williams et al., 2017), and soybean (Leung et al., 2020) within JULES-crop have been calibrated using observational data. Mathison et al. (2021) updated several rice and wheat parameters in JULES-crop, relying primarily on the literature, but did not account for O₃ effects. In this study, novel O₃-FACE experimental data were utilised to calibrate rice parameters in JULES-crop for the first time, improving its ability to assess O₃ impacts on rice growth.

JULES-crop utilises a flux-based approach to simulate the O₃ damage following Sitch et al. (2007). It assumes that the potential net photosynthesis A_p is suppressed by O₃:

$$A=A_{\mathrm{p}}F.$$

Here A is leaf-level net photosynthesis with the O₃ effects and F is the reduction factor:

$$F=\mathrm{1}-a\cdot max\left[F_{{\mathrm{O}}_{\mathrm{3}}}-F_{{\mathrm{O}}_{\mathrm{3}}\text{crit}},\mathrm{0}\right],$$

where $F_{{\mathrm{O}}_{\mathrm{3}}}$ represents instantaneous leaf uptake of O₃. $F_{{\mathrm{O}}_{\mathrm{3}}\text{crit}}$ and a are the plant-functional-type-specific threshold and sensitivity factor, respectively.

The O₃ flux $F_{{\mathrm{O}}_{\mathrm{3}}}$ ($\mathrm{nmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) is calculated as

$$F_{{\mathrm{O}}_{\mathrm{3}}}={\displaystyle \frac{\left[{\mathrm{O}}_{\mathrm{3}}\right]}{r_{\mathrm{a}}+\left[\frac{{\mathit{\kappa }}_{{\mathrm{O}}_{\mathrm{3}}}}{g_{\mathrm{l}}}\right]}}.$$

Here [O₃] (nmol m⁻³) is the molar O₃ concentration at the reference level, r_a (s m⁻¹) is the aerodynamic resistance and the boundary layer resistance between the leaf surface and reference level (Monin and Obukhov, 1954), ${\mathit{\kappa }}_{{\mathrm{O}}_{\mathrm{3}}}$ is the ratio of leaf conductance for O₃ to leaf conductance for water vapour (1.67), and g_l represents the leaf conductance for H₂O as a linear function of the photosynthetic rate (Cox et al., 1999):

$$g_{\mathrm{l}}=g_{\mathrm{l}}^{\ast }F,$$

where $g_{\mathrm{l}}^{\ast }$ is the leaf conductance in the absence of O₃ effects.

2.2 O₃-FACE experiments

The O₃-FACE experiment was conducted in Xiaoji, China ($\mathrm{32}\mathit{°}{\mathrm{35}}^{\prime }{\mathrm{5}}^{\prime \prime }$ N, $\mathrm{119}\mathit{°}{\mathrm{42}}^{\prime }{\mathrm{0}}^{\prime \prime }$ E), in 2012. It features four regular octagonal O₃-FACE fields (14 m in diameter) and four control fields, each covering an area of approximately 120 m². The experimental fields are spaced over 70 m apart to minimise the influence of O₃ release on neighbouring fields. Pipes positioned 50–60 cm above the crops released pure O₃ gas into each O₃-FACE field between 09:00 and 16:00 LT during the rice growing period. The mean daytime O₃ concentration during the experimental period was approximately 46 ppb under the elevated O₃ treatment, compared to 37 ppb in the ambient environment – an increase of around 25 %. The environmental conditions in the O₃-FACE and control fields were identical, except for the presence of O₃ pipes in the O₃-FACE fields. Samples from the O₃-FACE fields were collected from the field centre, at least 1.5 m away from the O₃ pipes, to ensure that the sampled rice had grown under stable O₃ conditions. Further details of the O₃-FACE system can be found in Wang et al. (2012).

The rice cultivar used was II You 084. The rice was planted on 30 May 2012 and reached maturity on 19 October 2012 in the ambient O₃ environment and 12 October 2012 in the elevated O₃ environment. During the growth period, key developmental stages, such as jointing and flowering, were recorded, and crop growth characteristics – including the dry biomass of leaves, stems, and panicles; LAI; and plant height – were measured at these stages to calibrate the model.

Three planting densities were employed during transplantation: low density (16 plants m⁻²), medium density (24 plants m⁻²), and high density (32 plants m⁻²). In addition to standard growth measurements, photosynthesis-related variables – including leaf temperature, internal leaf CO₂ concentration, stomatal conductance, and the photosynthesis rate (CO₂ assimilation rate) – were assessed using a LI-6400 portable photosynthesis system. After the rice reached maturity, 64 plants from each experimental field were harvested and dried to calculate the average rice yields.

2.3 FACE experiment for JULES-crop evaluation

Following calibration, observations of rice yields; height; and the dry weight of leaves, stems, and panicles from an independent FACE experiment were then used to evaluate the performance of JULES-crop. These additional field experiments were conducted in Danyang, China ($\mathrm{31}\mathit{°}{\mathrm{54}}^{\prime }{\mathrm{31}}^{\prime \prime }$ N, $\mathrm{119}\mathit{°}{\mathrm{28}}^{\prime }{\mathrm{21}}^{\prime \prime }$ E), and provided rice data for the 2022 and 2023 growing seasons. Two cultivars, Yangdao 6 and Wuyungeng 23, were transplanted on 20 July 2022 and 21 July 2023, respectively, and harvested between late October and early November. Yangdao 6 is an Indica rice cultivar, while Wuyungeng 23 belongs to the Japonica subspecies group, both of which represent the two major rice subspecies cultivated in China.

2.4 Data preparation

JULES-crop requires driving data, ancillary data, and control files to configure the model. Observations of hourly air pressure, specific humidity, air temperature, precipitation, wind speed, and shortwave radiation (SW) recorded during the O₃-FACE experiments were used as driving data. Diffuse radiation was calculated using a constant diffuse fraction in the model, with the default value of 0.4 applied in this study due to the absence of observational data. Surface downward longwave radiation (LW) was not measured in the O₃-FACE experiment and was instead estimated using an empirical model based on local observations (Chang and Zhang, 2019):

$$R_{\downarrow }=\mathit{\sigma }\cdot {\left(T_{\mathrm{a}}\right)}^{\mathrm{4}}\cdot \left[\text{clf}+\left(\mathrm{1}-\text{clf}\right)\cdot \left(a\cdot \ln \left({\displaystyle \frac{e_{\mathrm{a}}}{T_{\mathrm{a}}}}\right)+b\cdot \mathit{\phi }+c\right)\right].$$

Here R_↓ is the downward LW under all kinds of sky conditions (clear and cloudy), T_a is the air temperature; e_a is the water vapour pressure; φ is the relative humidity; σ is the Stefan–Boltzmann constant; a, b, and c are the empirical coefficients (Table 1); and clf is the cloud modification factor, set to 0 under clear-sky conditions:

$$\text{clf}=\mathrm{1}-K_{\mathrm{t}}.$$

Here K_t is the clearness index, which was calculated as follows:

$$K_{\mathrm{t}}={\displaystyle \frac{H_{\mathrm{m}}}{H_{\mathrm{0}}}},$$

where H_m represents the hourly measured solar radiation and H₀ denotes the hourly extraterrestrial solar radiation. Detailed calculation for H₀ can be found in Kumar and Umanand (2005).

Table 1 Empirical coefficients used in the longwave radiation model.

For the ancillary data, soil property values were extracted from the ancillary dataset used in the HadGEM2-ES model, which also underpins global simulations (Osborne et al., 2015). Another crucial factor influencing crop growth, the annual average CO₂ concentration, was set based on data provided by the Global Monitoring Laboratory (GML) of the National Oceanic and Atmospheric Administration (NOAA).

The weather station for the evaluation experiments provided only daily temperature and precipitation data. Consequently, additional meteorological variables, including wind, humidity, and LW, were sourced from the ECMWF Reanalysis v5 (ERA5) dataset. However, the ERA5-generated SW for 2022 and 2023 disrupted the JULES-crop simulations leading to unrealistically high leaf area index (LAI) values (exceeding 15). The overestimation of SW in ERA5 has been widely reported, with studies attributing it to the omission of aerosol variations and a limited capacity to simulate clouds and water vapour, resulting in an overestimation of hourly SW in China by approximately 73.95 W m⁻² (He et al., 2021; Jiang et al., 2020; Tong et al., 2023; Li et al., 2023). To address this, SW was bias-corrected using observations from the O₃-FACE experiment conducted in 2012.

Additionally, O₃ concentration observations were unavailable for the evaluation experiments. Hourly O₃ data from the nearest station of the China National Environmental Monitoring Centre (https://www.cnemc.cn/, last access: 20 December 2024) were used instead. Aside from these driving data, e.g. weather variables, O₃ concentrations, CO₂ concentrations, and crop stage dates, the evaluation simulations applied the same settings and parameters as those used in the calibration.

3 Results

3.1 Calibration

All parameters calibrated using the O₃-FACE experiment are listed in Tables 2 and 3. The calibration process for rice involved four main steps. First, leaf-level simulations were calibrated by fitting simulated photosynthesis rates with observed values. The nitrogen content in leaves, stems, and roots was obtained from observations and the literature. Observed leaf temperature, internal CO₂ concentration, and stomatal conductance were used as model inputs. Photosynthesis-related parameters were adjusted based on discrepancies between observed and simulated photosynthesis rates. Notably, O₃ damage was not considered during this step.

Table 2 Calibrated plant functional type (PFT) parameters representing rice.

Table 3 Calibrated crop-related parameters representing rice.

* These parameters were spatially varying in Osborne et al. (2015).

Second, canopy-level simulations were calibrated by determining the rice growth rate and partitioning of assimilated carbon. Air temperature data were used to calculate the accumulated temperature required for rice growth stages, and the allocation of carbon to various carbon pools was also defined during this phase.

Third, model simulations were evaluated against observed LAI and crop height following the calibration of crop physiology parameters. Lastly, rice yields were compared with observations under both ambient and elevated O₃ concentrations.

The calibration process involved iteratively adjusting parameters manually until the model simulations fell within the range of observed values. Additional adjustments were made to refine results, aiming to align them more closely with the central tendency of the observations. Although the number of simulations was constrained by computational limitations, the process successfully achieved agreement with all available observations, ensuring no discrepancies remained. While finer and finer incremental adjustments were not feasible due to computational limitations, the approach effectively balanced precision and generalisation, capturing the essential crop observations without overfitting.

3.1.1 Photosynthesis

The potential leaf-level photosynthesis, unaffected by water stress and O₃ effects, is calculated based on three potentially limiting rates: the Rubisco-limited rate (W_c), the light-limited rate (W_l), and the rate of transport of photosynthetic products (W_e) for C₃ plants, as detailed in Clark et al. (2011).

Figure 2 Leaf nitrogen concentration (kg N per kg C) (a) and the ratio of root nitrogen concentration to leaf nitrogen concentration (b). The dashed grey line represents the values selected for the simulation, while the dots indicate the observed values.

Following Farquhar et al. (1980) and Collatz et al. (1991), several parameters in the photosynthesis scheme are temperature-dependent, including the maximum rate of Rubisco carboxylation, V_m ($\mathrm{mol}{\mathrm{CO}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$), which is critical for calculating both W_c and W_e. V_cmax is calculated assuming an optimal temperature range defined by T_upp and T_low.

$$V_{\text{cmax}}={\displaystyle \frac{V_{\text{cmax25}}{f}_T\left(T_{\mathrm{c}}\right)}{\left[\mathrm{1}+e^{\left\{\mathrm{0.3}\left(T_{\mathrm{c}}-T_{\text{upp}}\right)\right\}}\right]\left[\mathrm{1}+e^{\left\{\mathrm{0.3}\left(T_{\text{low}}-T_{\mathrm{c}}\right)\right\}}\right]}},$$

where V_cmax25 represents the maximum rate of carboxylation of the enzyme Rubisco at 25 °C and is assumed to be linearly dependent on the leaf nitrogen concentration. For the C₃ crop, V_cmax25=n_en_l, where n_e is the scale factor and n_l is the leaf nitrogen concentration (kg N per kg C). T_c is the leaf temperature in °C, T_upp and T_low are PFT-dependent parameters, and f_T depends on the parameter q_10,leaf, the factor by which plant respiration increases by a 10 °C increase in temperature:

$$f_T=q_{\mathrm{10},\text{leaf}}^{\mathrm{0.1}(T_{\mathrm{c}}-\mathrm{25})}.$$

Changes in PFT parameters primarily influence the simulations of the photosynthesis rate, which in turn affects the accumulation of carbon in rice. In JULES-crop, the photosynthesis process is closely linked to the nitrogen content of the crop. Leaf nitrogen concentration (n_l) is a key factor impacting the photosynthesis rate and was estimated based on literature sources (Fig. 2a). As leaf nitrogen concentration declines from the vegetative to the ripening stage, the rice plant's capacity for carbon accumulation diminishes.

The ratio of the nitrogen content of roots relative to leaves (μ_rl) was also derived from the literature (Fig. 2b). This ratio determines the nitrogen content in the roots, which further influences the respiration rate. The maturity stage was excluded when calculating the average values for each stage. The values presented in Fig. 2 were collected from peer-reviewed studies conducted across China over the past 20 years (listed in the Supplement), encompassing several rice cultivars grown in major rice-producing regions.

The ratio of the nitrogen content of stems to that of leaves (μ_sl) was determined from the O₃-FACE observations. The ratio varied across growth stages, reaching its highest value during the maturity stage (Fig. 3). This is because at maturity the leaves consist solely of yellow leaves, which have lower nitrogen content compared to the green leaves present during earlier stages. The calibrated μ_sl is the average value during the tillering, jointing, and heading stages.

Figure 3 Ratio of the stem nitrogen concentration to the leaf nitrogen concentration. The dashed grey line and the dots show the values for the simulation and observations, respectively.

The simulations of the net leaf photosynthesis rate, using the default parameters from Osborne et al. (2015), underestimated the observed values (Fig. 4a). Several parameters, including n_l, n_e, f_dr, T_upp, and q_10,leaf, were calibrated to make the simulation results agree better with observations. The standard photosynthesis model assumes that the upper temperature limit for C₃ crops is 36 °C. However, when the temperature exceeded 36 °C, the simulated photosynthesis rates were still underestimated (Fig. 4b). This suggests that temperatures above 36 °C should be increased to 38 °C to obtain improved agreement with observations, as shown in Figs. 4c and 5.

Figure 4 Simulated photosynthesis rate ($\mathrm{\mu }\mathrm{mol}{\mathrm{CO}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) using parameters before calibration (a) and after calibration without (b) or with (c) changing the upper temperature limitation parameter (T_upp). The dashed line is the 1:1 line.

Figure 5 The coloured lines are simulated temperature responses of the photosynthesis rate using the mean value of the observed intercellular CO₂ concentration of leaves (Ci) and calibrated (38 °C) or default (36 °C) T_upp. The filled dots and open circles represent the observations used in this study and simulations generated by calibrated parameters, respectively. The error bars are taken from five independent studies (Table S1 in the Supplement) that span multiple rice cultivars, nitrogen regimes, CO₂ levels, and light intensity levels.

Figure 5 shows that the simulated leaf photosynthetic rate starts to decrease at approximately 30 °C using the calibrated temperature parameters, while the simulated curves using the default T_upp from Osborne et al. (2015) reached the optimum temperature at about 29 °C. The exact optimum temperature for simulations varied with the intercellular CO₂ concentration of leaves (Ci). According to the experimental data collected from the literature, the optimum temperature should be around 30 °C, depending on the environmental conditions such as the nitrogen content of leaves, light intensity, and CO₂ concentration as well as growth stages. After calibration, the response of the leaf photosynthetic rate to leaf temperature was closer to observations from both this study and the literature.

3.1.2 Rice development and assimilate partitioning

The development status of rice is closely linked to its phenological progression and is represented by the development index (DVI). The DVI increases as the ratio of accumulated thermal time to the prescribed thermal time for each developmental phase rises. Initially, the DVI is set to −1 at sowing, increases to 0 at emergence, completes accumulation before flowering at a value of 1, and reaches a value of 2 at maturity.

Once rice is sown, its developmental rate, defined by the DVI, depends on the prescribed thermal time, which includes the thermal time between sowing, emergence, flowering, and maturity stages (Osborne et al., 2015). The thermal time (T_eff) can be calculated as follows:

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

where T, T_b, T_o, and T_m are air temperature, base temperature (8 °C), optimum temperature (30 °C), and maximum temperature (42 °C), respectively, which are the values from Osborne et al. (2015).

The changes in the value of DVI during the simulation are determined by

$${\displaystyle \frac{\mathrm{d}\text{DVI}}{\mathrm{d}t}}=\left\{\begin{array}{ll}{T}_{\text{eff}}/T_{\text{emr}},& \text{for }-\mathrm{1}\le \text{DVI}<\mathrm{0}\\ T_{\text{eff}}/T_{\text{veg}},& \text{for }\mathrm{0}\le \text{DVI}<\mathrm{1}\\ T_{\text{eff}}/T_{\text{rep}},& \text{for }\mathrm{1}\le \text{DVI}<\mathrm{1}\end{array}\right.,$$

where T_emr, T_veg, and T_rep represent the thermal time intervals between sowing and emergence, emergence and flowering, and flowering and maturity, respectively.

The field experiment recorded the dates for sowing, transplanting, panicle initiation, heading, and maturity when collecting samples. In China, most rice is grown in puddled fields after transplanting (Wang et al., 2017). Before transplanting, rice is cultivated in nurseries and is not moved to the field until it has developed five or six leaves. The prescribed thermal time was estimated based on the calculated thermal time from the observations (Table 4). The observed development stages and crop characteristics are used to determine the thermal time required for the model, ensuring that the following conditions are met: the model's predicted maturity stage coincides with the actual timing observed in the experiment and the DVI of crop characteristics from simulations agrees with the observations. For example, the transplanting stage falls within the vegetative phase, so the DVI of observations should fall within the range of 0 to 1.

Table 4 Thermal time of rice between transplanting and maturity.

Once the development rate was determined, the accumulated net primary productivity (NPP) of each time step was partitioned into four main carbon pools: root, stem (including structural stem and stem reserves), leaves, and the harvest pool (including yellow leaves and harvested organs, which are panicles for rice).

The partition coefficients (p) are calculated as follows (Osborne et al., 2015):

$$\begin{array}{rl}{\displaystyle }{p}_{\text{root}}& {\displaystyle }={\displaystyle \frac{e^{-{\mathit{\alpha }}_{\text{root}}+{\mathit{\beta }}_{\text{root}}\text{DVI}}}{e^{-{\mathit{\alpha }}_{\text{root}}+{\mathit{\beta }}_{\text{root}}\text{DVI}}+e^{-{\mathit{\alpha }}_{\text{stem}}+{\mathit{\beta }}_{\text{stem}}\text{DVI}}+e^{-{\mathit{\alpha }}_{\text{leaf}}+{\mathit{\beta }}_{\text{leaf}}\text{DVI}}+\mathrm{1}}},\\ {\displaystyle }{p}_{\text{stem}}& {\displaystyle }={\displaystyle \frac{e^{-{\mathit{\alpha }}_{\text{stem}}+{\mathit{\beta }}_{\text{stem}}\text{DVI}}}{e^{-{\mathit{\alpha }}_{\text{root}}+{\mathit{\beta }}_{\text{root}}\text{DVI}}+e^{-{\mathit{\alpha }}_{\text{stem}}+{\mathit{\beta }}_{\text{stem}}\text{DVI}}+e^{-{\mathit{\alpha }}_{\text{leaf}}+{\mathit{\beta }}_{\text{leaf}}\text{DVI}}+\mathrm{1}}},\\ {\displaystyle }{p}_{\text{leaf}}& {\displaystyle }={\displaystyle \frac{e^{-{\mathit{\alpha }}_{\text{leaf}}+{\mathit{\beta }}_{\text{leaf}}\text{DVI}}}{e^{-{\mathit{\alpha }}_{\text{root}}+{\mathit{\beta }}_{\text{root}}\text{DVI}}+e^{-{\mathit{\alpha }}_{\text{stem}}+{\mathit{\beta }}_{\text{stem}}\text{DVI}}+e^{-{\mathit{\alpha }}_{\text{leaf}}+{\mathit{\beta }}_{\text{leaf}}\text{DVI}}+\mathrm{1}}},\\ {\displaystyle }{p}_{\text{harv}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{e^{-{\mathit{\alpha }}_{\text{root}}+{\mathit{\beta }}_{\text{root}}\text{DVI}}+e^{-{\mathit{\alpha }}_{\text{stem}}+{\mathit{\beta }}_{\text{stem}}\text{DVI}}+e^{-{\mathit{\alpha }}_{\text{leaf}}+{\mathit{\beta }}_{\text{leaf}}\text{DVI}}+\mathrm{1}}}.\end{array}$$

Six parameters, α_root, α_stem, α_leaf, β_root, β_stem, and β_leaf, determine the partitioning process during the whole growth period. And the NPP accumulated through photosynthesis for each time step is distributed to the four carbon pools according to the partition coefficients.

The parameters for carbon distribution were calibrated based on the dry weights of the stem, leaves, and panicles from field experiments. Since root carbon is not included in the observations, its partitioning value is estimated as a fraction of rice yield. Liu et al. (2023) suggested that the ratio of root dry weight to grain yield is approximately 0.13, although it can vary depending on the cultivar and nitrogen application rate. Figure 6 shows the fraction of accumulated NPP partitioned into the different carbon pools using the calibrated parameters.

Figure 6 Fraction of daily accumulated net primary productivity partitioned into roots (blue), stems (orange), leaves (green), and harvested parts (red) of the crop as a function of the development index (DVI; 0 = emergence, 1 = flowering, 2 = maturity) for rice. The dotted black line is the fraction based on parameters used in Osborne et al. (2015).

Figure 7 Specific leaf area against the development index. Coloured symbols indicate observations, and the colour shows the data from different experiment fields. The dotted black line and the solid black line show the fit using parameters from Osborne et al. (2015) and our calibrated parameters, respectively. Note that various symbols correspond to successive sampling dates from the same experimental field, thereby illustrating the temporal progression of the observations.

Figure 8 Stem dry weight against crop height. Coloured symbols are observations, and the colour shows the data from O₃-FACE experiment. The dashed black line and the solid black line show the fit using parameters from Osborne et al. (2015) and the calibrated parameters, respectively. Note that various symbols correspond to successive sampling dates from the same experimental field, thereby illustrating the temporal progression of the observations.

The accumulated carbon in different carbon pools directly affects the biomass of various rice organs. The model calculates carbon accumulation and distribution, so the fractions of carbon to dry matter in the root, stem, and leaf (f_C, root, f_C, stem, and f_C, leaf) must be defined prior to running the model. The values used in our calibrated simulations were taken from the observations and are listed in Table 2, along with the default values from Osborne et al. (2015). The value of the carbon fraction impacts the root growth, crop height, and LAI.

Figure 9 Carbon-to-dry-matter ratio of the leaf, stem, root, and panicle during different crop development stages, where the average means the value collected from the literature, which only provided an average value for all stages during the rice growth. The dashed green, yellow, and blue lines represent the value prescribed in the model for the fraction of carbon to dry matter in the root, stem, and leaf, respectively.

Figure 10 Leaf, stem, and total aboveground carbon against the day of year under ambient and elevated O₃ conditions. Box plots are observations, whereas lines show the simulation results using parameters from Osborne et al. (2015) (grey) and calibrated (green) parameters under ambient O₃ conditions, including high (blue) and low (orange) O₃ sensitivity under elevated O₃ conditions, with units of g m⁻².

Table 5 O₃ parameters calibrated for high and low sensitivity to O₃ damage.

3.1.3 LAI and crop height

LAI is an important attribute of crops, reflecting their capacity for carbon accumulation. In JULES-crop, LAI is linked to the leaf carbon pool (Osborne et al., 2015):

$$\text{LAI}={\displaystyle \frac{C_{\text{leaf}}}{f_{\text{C,\hspace{0.17em}leaf}}}}\text{SLA}.$$

Here C_leaf indicates the amount of carbon in leaves, f_C, leaf represents the carbon fraction of dry matter in leaves, and SLA is the specific leaf area (m⁻² leaf kg⁻¹):

$$\text{SLA}=\mathit{\gamma }(\text{DVI}+\mathrm{0.06}{)}^{\mathit{\delta }},$$

where γ and δ are determined by fitting the curve between DVI and SLA (De Vries et al., 1989) from observations (Fig. 7).

As green leaves begin to turn yellow, leaf senescence starts and is represented by the parameter DVI_sen. The change from green to yellow signals the transition of carbon from the leaf carbon pool to the harvest carbon pool. The transition rate is simulated by reducing C_leaf by a specific fraction (De Vries et al., 1989):

$$C_{\text{harv}}=C_{\text{harv}}+\mathit{\mu }{\left(\text{DVI}-{\text{DVI}}_{\text{sen}}\right)}^{\mathit{\nu }}\cdot C_{\text{leaf}},$$

where μ and ν were determined by fitting the declining trend of carbon in green leaves following leaf senescence. The simulation results are presented in Sect. 3.1.4.

The calculation of crop height (h) depends on the amount of carbon in the stem (C_stem) (Hunt, 2012):

$$h=\mathit{\kappa }{\left({\displaystyle \frac{C_{\text{stem}}}{f_{\text{C,\hspace{0.17em}stem}}}}\right)}^{\mathit{\lambda }},$$

where f_C, stem represents the carbon fraction of dry matter in the stem and κ and λ are determined by fitting the relationship between h and the stem dry matter of stems, which is equal to $\frac{C_{\text{stem}}}{f_{\text{C,\hspace{0.17em}stem}}}$ (Fig. 8).

Similarly to leaf senescence, the carbon stored in the stem reserves is mobilised into the harvest carbon pool at a rate of 10 % d⁻¹ once the partition coefficient for stems drops below 0.01 (De Vries et al., 1989).

$$C_{\text{harv}}=C_{\text{harv}}+\mathrm{0.1}\cdot \mathit{\tau }{C}_{\text{stem}},$$

where τ represents the fraction of stem growth partitioned to reserve carbon.

The observations did not include components of the carbon fraction, such as C_leaf and C_stem, required for the model simulation; therefore, these values were sourced from peer-reviewed rice field studies (listed in the Supplement). Some studies evaluated varied stressors or environmental treatments. Thus, to ensure consistency with calibration, only the control-plot values under default (unstressed) conditions were used. All the literature data were derived from rice field experiments conducted in China, involving several rice cultivars to enhance representativeness (Fig. 9). The carbon content of panicles was also obtained from the literature and combined with the carbon in yellow leaves during the ripening phase to calculate the total carbon in the harvest pool. Additionally, the fractions of carbon to dry matter were used to compare the simulation results with the observations, which only provided dry biomass data for rice.

3.1.4 Comparison with O₃-FACE experiments

Figure 10 illustrates the changes in the main carbon pools throughout the entire growing period. The accumulated carbon was reduced under elevated O₃ conditions, highlighting the detrimental impact of O₃ on crop growth. At the maturity stage, total aboveground carbon under elevated O₃ was 22 %–29 % lower compared to ambient O₃ conditions, as shown in the observations (Fig. 10e and f). Carbon levels in both the leaf and stem exhibited a similar decreasing trend due to the O₃-induced damage to the photosynthesis process and carbon accumulation. The simulations closely matched the observations, using the average carbon-to-dry-matter fraction for different growth stages to convert observed data into carbon weights (Fig. 9). It is important to note that the carbon fraction varies with cultivar and growing environment. To align the model results, which are based on carbon weight instead of dry weight, with the observed data, the average carbon-to-dry-matter ratio across all stages was applied.

There are two parameters in the simulations that directly relate to the impact of O₃ on the rice (Clark et al., 2011; Sitch et al., 2007) (Table 5). The reduction in the net photosynthesis rate was determined by the value of the instantaneous leaf uptake of O₃ above the threshold $F_{{\mathrm{O}}_{\mathrm{3}}\text{crit}}$, multiplied by a sensitivity parameter a (Pleijel et al., 2004). Observations using three planting densities of rice observations were used to calibrate the model. As can be seen from Fig. 11, the high sensitivity and low sensitivities coincided with the upper and lower boundaries of relative yield (RY), which is calculated as follows:

$$\text{RY}={\displaystyle \frac{Y_{{\mathrm{O}}_{\mathrm{3}}}}{Y_{\mathrm{0}}}},$$

where $Y_{{\mathrm{O}}_{\mathrm{3}}}$ represents the crop yield including O₃ damage and Y₀ represents the crop yield with no effects of O₃.

Figure 11 Relative yield against AOT40 (ppm h). Coloured lines show the relative yield of rice planted in high (green), medium (orange), and low (blue) density. The grey lines show the simulations of relative yield with high (dotted) and low (dashed) sensitivity to O₃ damage.

In Fig. 11, AOT40 was used to represent the O₃ concentrations in the environment:

$$\text{AOT}\mathrm{40}=\sum _{i=\mathrm{1}}^n\left({\left[{\mathrm{O}}_{\mathrm{3}}\right]}_i-\mathrm{0.04}\right),$$

where [O₃]_i stands for the hourly O₃ concentration level (unit: ppm h) during daylight hours (08:00–19:59) and n represents the total hours of the growing season.

Figure 12 illustrates the height and LAI of rice under both elevated and ambient O₃ conditions. The difference in LAI and height between these two environments underscores the negative impact of O₃ on rice carbon accumulation. Post-calibration, the simulations for both LAI and height align well with observational data (Fig. 11a and c). Prior to the new calibration, simulations with default parameters from Osborne et al. (2015) significantly underestimated both LAI and height largely due to the underestimated photosynthesis rate (Fig. 4). This underestimation led to reduced carbon assimilation and storage, resulting in insufficient carbon allocation to stems and leaves, which directly impacted LAI and height. It is worth noting that all plots comparing simulation and observation begin after the model's initialisation phase.

Figure 12 Crop height (cm) and green leaf area index (LAI) are shown versus the day of year under ambient and elevated O₃ conditions. Box plots show observations, whereas lines show the simulation results using parameters from Osborne et al. (2015) (grey) and calibrated (green) parameters under ambient O₃ conditions, including high (blue) and low (orange) O₃ sensitivity under elevated O₃ conditions.

3.2 Evaluation

Figure 13 compares simulated and observed values of leaf carbon, stem carbon, total aboveground biomass, and rice height for the years 2022 and 2023, based on data from an independent FACE experiment (see Sect. 2.2). The observations were limited to heading and maturity stages. These observations were compared to our newly calibrated JULES-crop model simulations using these FACE observations. The calibrated O₃-damage parameters were applied to model the impact of O₃ on rice biomass and carbon content.

Figure 13 Leaf, stem, total aboveground biomass, and crop height against the day of year for 2022 and 2023. Box plots are observations, and coloured lines show the simulation results using low (orange) and high (blue) O₃ sensitivity, respectively.

The simulated stem carbon was marginally lower than the average observed values (Fig. 13c and d), while total aboveground biomass was overestimated when using low O₃ sensitivity parameters (Fig. 13e and f). These variations can be attributed to differences in the carbon allocation between the calibration and evaluation experiments. The seeding depth notably influenced stem weight since stems thickened nearer the root, and only aboveground stems were harvested and measured. Consequently, deeper seeding resulted in a smaller fraction of stem biomass relative to total aboveground biomass (Gong et al., 2023). This slight underestimation of stem carbon was also evident in the simulation of crop height, which was similarly affected by seeding depth.

The total biomass observed in the evaluation experiment surpassed that measured in the O₃-FACE experiment, particularly with a notably larger stem weight. Crop parameters were calibrated using data from the O₃-FACE experiment, but differences in agronomic practices across experiments may have introduced uncertainties.

While the simulated crop height fell within the range of observed values, it was marginally lower than the average measured height (Fig. 13g and h). Despite variations in seeding practices between the calibration and evaluation field experiments, the carbon distribution and levels aligned well with the observations. Overall, JULES-crop demonstrated the ability to accurately predict rice growth and carbon allocation across various carbon pools.

3.3 Limitations

While this study provides a rice model calibration based on the novel O₃-FACE experiments, several limitations must be acknowledged. The calibrated O₃ parameters influence modelled net photosynthesis, biomass, and yield through the control on stomatal uptake and instantaneous photosynthesis. Due to limited O₃-FACE observations, our calibration did not represent differences in O₃ sensitivity due to rice cultivar.

Additionally, the calibrated thermal time was specific to a particular location and should be recalculated using local air temperature and rice phenology data if simulations are performed for other regions. For example, the evaluation experiment conducted in 2023 in a nearby county exhibited a relatively high thermal time compared to the calibration experiment, primarily due to the longer growth duration. Rice growth in the 2022 and 2023 evaluation experiments was severely affected by crop pests and diseases at the maturity stage, leading to significant yield loss. As a result, only crop growth characteristics were used to validate the model.

Furthermore, although the model was calibrated and evaluated using independent experimental data, directly applying the parameters to global simulations may introduce significant uncertainties. As such, global simulations using the parameters derived in this study should incorporate further evaluations to verify model performance (Müller et al., 2017).

4 Conclusions

This study marks a significant advancement in modelling rice growth and O₃ effects by providing the first calibration of the JULES-crop model using rice-specific data from O₃-FACE experiments. These experiments offer a realistic field setting to assess the impacts of O₃ on crops, addressing the limitations of alternative setups such as OTCs by simulating more natural environmental conditions. Initial simulations with the default rice parameters in JULES-crop revealed substantial underestimation of carbon accumulation throughout the growth cycle. Calibration using the most recent O₃-FACE data significantly improved the model's ability to replicate rice physiology, phenology, yield, and O₃ sensitivity.

The calibration process involved adjusting key parameters to align simulations with observed data, including leaf area indices; crop height; yield; and the biomass of leaves, stems, and panicles. The model was refined to accurately represent yield reductions caused by elevated O₃ levels. Evaluation against independent field experiments demonstrated good agreement between simulated outcomes and observed results, affirming the model's robustness.

This study deepens our understanding of O₃'s impact on rice production and delivers a newly calibrated model suitable for assessing future climate scenarios and O₃ effects. The study lays the groundwork for future agricultural research aimed at mitigating O₃-induced yield losses, providing a valuable framework for enhancing food security as O₃ levels continue to rise.