O₃–precursor relationship over multiple patterns of timescale: a case study in Zibo, Shandong Province, China

In this study, we developed an approach that integrated multiple patterns of timescale for box modeling (MCMv3.3.1) to better understand the O₃–precursor relationship at multiple sites and through continuous observations. A 5-month field campaign was conducted in the summer of 2019 to investigate the ozone formation chemistry at three sites in a major prefecture-level city (Zibo) in Shandong Province of northern China. It was found that the relative incremental reactivity (RIR) of major precursor groups (e.g., anthropogenic volatile organic compounds (AVOCs), NO_x) was overall consistent in terms of timescales changed from wider to narrower (four patterns: 5-month, monthly, weekly, and daily) at each site, though the magnitudes of RIR varied at different sites. The time series of the photochemical regime (using RIR${}_{{\mathrm{NO}}_x}$ $/$ RIR_AVOC as an indicator) in weekly or daily patterns further showed a synchronous temporal trend among the three sites, while the magnitude of RIR${}_{{\mathrm{NO}}_x}$ $/$ RIR_AVOC was site-to-site dependent. The derived RIR ranking (top 10) of individual AVOC species showed consistency between three patterns (i.e., 5-month, monthly, and weekly). It was further found that the campaign-averaging photochemical regimes showed overall consistency in the sign but non-negligible variability among the four patterns of timescale, which was mainly due to the embedded uncertainty in the model input dataset when averaging individual daily patterns into different timescales. This implies that utilizing narrower timescales (i.e., daily pattern) is useful for deriving reliable and robust O₃–precursor relationships. Our results highlight the importance of quantifying the impact of different timescales to constrain the photochemical regime, which can formulate more accurate policy-relevant guidance for O₃ pollution control.

1 Introduction

Since 2013, the ambient PM_2.5 concentration in China has dramatically declined following the implementation of the Clean Air Action (Lu et al., 2018; Y. Wang et al., 2020; Zhang et al., 2019). However, national ground surface ozone concentrations increased over the same period (Xue et al., 2020) and became a major air quality problem that needed to be addressed in China (Li et al., 2019; Wang et al., 2019). It is well-known that ground surface ozone is formed mainly by complex nonlinear photochemical oxidation of volatile organic compounds (VOCs) in the presence of nitrogen oxides (NO_x= NO + NO₂) and sunlight (Blanchard, 2000; Hidy, 2000; Kleinman, 2000), which adversely influences human health, vegetation, and crops (Brunekreef and Holgate, 2002; Vingarzan, 2004).

Given the complex non-linear relationship between O₃ formation and its precursors (VOCs and NO_x), challenges in mitigating its severity lie primarily in comprehensively understanding the O₃–precursor relationship (Su et al., 2018a; Tan et al., 2018a). It is commonly recognized that regional-scale air quality models and the 0-D box model are two mainstream approaches to investigate the increasingly severe ozone problem (Blanchard, 2000; Cardelino and Chameides, 1995; Hidy, 2000; Liu et al., 2019). Unlike the complicated 3-D air quality models, the 0-D box model is an observation-based model that is implemented with a gas-phase chemical mechanism; it has been widely used to diagnose O₃–precursor relationships in various locations (Liu et al., 2021a; Sun et al., 2016; Tan et al., 2019; Xue et al., 2014a; Yu et al., 2020a). Some previous studies (Li et al., 2021; Lu et al., 2010a; Sicard et al., 2020; Yu et al., 2020b) have reported large variability in O₃–precursor relationships on spatiotemporal scales in many cities of China, which indicates great challenges for current O₃ pollution control (Y. Wang et al., 2017; Xue et al., 2014b).

Table 1 summarizes the published studies of O₃–precursor relationships using the 0-D box model (implemented with different gas-phase chemical mechanisms) with diversified patterns of timescale in many places in China. The observational period in most previous studies was short term (i.e., less than one month), while medium-term (i.e., from one to several months), and long-term (i.e., multiple years) periods were limited. As shown in Table 1, we find that model input datasets with different timescales have been employed in previous studies to identify the campaign-averaging O₃ formation regime, but there is a lack of comparison among these different timescales. We also find that more than half of the studies use the averaged diurnal patterns as box model input, which is particularly common for medium- and long-term measurements. For example, a 10-year long-term observational study by Y. Wang et al. (2017) adopted a monthly pattern of timescale for model simulation for the reason of saving computing resources; it also revealed a substantial temporal variability in the O₃–precursor relationship. In addition, it is believed that long-term (measurements of at least several months) and multiple-site continuous online measurements can provide an opportunity to develop O₃ control strategies more comprehensively over a wider spatiotemporal scale (Li et al., 2021; Y. Wang et al., 2017; T. Wang et al., 2017). However, such measurements have been quite rare in China, limiting the present understanding of O₃–precursor relationships (Lu et al., 2019; T. Wang et al., 2017).

Table 1 Summary of relevant published 0-D box model studies in China.

^a Number of days for modeling the patterns of timescale denotes that which was simulated by the box model. ^b Peking University. ^c Shenzhen Yanjiusheng Yuan. ^d Hubei Provincial Environmental Monitoring Center. ^e Nanjing University of Information Science & Technology. ^f Environmental Monitoring Building. ^g Environmental Protection Bureau.

In this study, a 5-month field campaign was conducted in the summer of 2019 to investigate the ozone formation chemistry at three sites in Zibo, a major prefecture-level Chinese city in Shandong Province. According to our measurements at the three sites in Zibo, the averaged O₃ concentration during the whole observational period was around 50 ppbv, while the daily maximum of O₃ concentrations for some extremely polluted periods were nearly 120–150 ppbv (see details in Sect. 3.1). Here, we developed an approach that integrated multiple patterns of timescale for box model simulation, which aimed to illustrate the non-linearity of O₃–precursor relationships driven by their actual daily, weekly, and monthly variability. Our results can be conducive to interpreting variations of O₃–precursor relationships over a wider spatiotemporal scale, and they provide implications for developing more precise and constrained O₃ control strategies in other regions.

2 Methods

2.1 Study sites and measurements

Field measurements were conducted in a major prefecture-level city (Zibo), which is in the middle of Shandong Province, northern China, from 1 May to 30 September 2019. Figure S1 in the Supplement shows the surrounding environment and geographical locations at the three sampling sites; a detailed description of the Tianzhen (TZ), Beijiao (BJ), and Xindian (XD) sites can be found in our previous study (Li et al., 2021). Briefly, TZ contains a mixture of crude oil processing and operation stations and farming areas and is classified as a suburban area; XD contains a mixture of residential and heavy industrial zones and is considered to be a suburban area; BJ is in the urban area of Zibo.

Typical inorganic gases of O₃, NO, NO₂, CO, and SO₂ were measured using online commercial gas analyzers (Thermo Scientific 49i, 42i, 48i, and 43i, USA) at the three sites. Following the Chinese meteorological monitoring regulation (GB/T 35221-2017), we continuously monitored the meteorological parameters (i.e., temperature, relative humidity, UV-A solar radiation, precipitation, wind speed, and wind direction) at the three sites (Li et al., 2021). Two online GC systems (gas chromatography–flame ionization detector, GC–FID, Thermo Scientific GC5900) were deployed at TZ and BJ respectively to measure VOC species. For C₂-C₅ VOCs, desorption and separation were performed using a GC with pre-concentration on a combination of two columns, followed by an FID detector. For C₆-C₁₂ VOCs, the air sample was pre-concentrated on Tenax GR cartridges and subsequently separated by chromatographic column, then it was detected by another FID detector. Similarly, one online system (gas chromatography–flame ionization detector–photoionization detector, GC–FID–PID, Syntech Spectras GC 955-615/815) was deployed at the XD site. For C₂-C₆ VOCs, the hydrocarbons were concentrated on a Tenax GR carrier then thermally desorbed and separated on a DB-1 column before finally being detected by the FID and PID. For C₆-C₁₂ VOCs, the air sample was concentrated on a Carbosieve SIII carrier at 5^∘ then thermally desorbed and separated on a combination of two columns; the FID and PID detectors were employed for subsequent detection. These systems measured 55 VOC species at a 1 h resolution; more detailed descriptions can be found elsewhere (Chien, 2007; Jiang et al., 2018; Xie et al., 2008).

Table S1 summarized the limit of detection, accuracy, and precision of the instruments at the three sites, and all the measurement instruments were regularly subjected to service checking and maintenance during the whole campaign. Unfortunately, we did not conduct the inter-comparison between the GC–FID and GC–FID–PID instruments at the same site due to practical reasons, as these VOC instruments were separately deployed at the three different sites for continuous routine operation. To ensure the quality assurance and quantity control (QA–QC) of online VOC measurement, two five-point calibrations (i.e., 2, 4, 6, 8, 10 ppbv, dilution from one cylinder) for standard gases with 55 VOC species (Linde Co., Ltd, USA) were carried out in May and August of 2019 at the three sites. Table S2 showed that the calibration linearity (R²) of all measured VOCs was nearly 0.9990. Additionally, a single-point calibration (i.e., 6 ppbv) was regularly performed every month during the whole campaign. As shown in Fig. S2 (a case from TZ), the retention time, peak fitting, and baseline of the chromatogram were manually checked and adjusted on a daily basis.

2.2 0-D box model and design of four patterns of timescale

The 0-D box model integrated with the latest Master Chemical Mechanism MCMv3.3.1 (http://mcm.york.ac.uk, last access: 27 January 2023) has been widely utilized in many regions (He et al., 2019; Jenkin et al., 2015; Liu et al., 2019; Whalley et al., 2021). Unlike the lumped chemical mechanisms such as CB05 (Y. Wang et al., 2017; Yarwood et al., 2005), CB6 (Yarwood et al., 2010), RACM and RACM2 (Goliff et al., 2013; Stockwell et al., 1997, 2020), and SAPRC-07 (Carter, 2010), the MCMv3.3.1 is a near-explicit chemical mechanism consisting of over 5800 species and 17 000 reactions (Jenkin et al., 2015; Saunders et al., 2003), which can be used to describe the gas-phase chemistry (i.e., in situ photochemistry). In this study, the box model (based on the Framework for 0-D Atmospheric Modeling, F0AM; Wolfe et al., 2016) was applied and constrained by the mean diurnal profiles of meteorological data (i.e., temperature, relative humidity, and photolysis rates), 4 inorganic gases (i.e., SO₂, CO, NO, and NO₂), and 45 speciated VOCs (in MCMv3.3.1 species list; see Table S3). Since measured photolysis rates (J values) were not available, the measured UV-A solar radiation was used to scale the photolysis rates calculated from the Tropospheric Ultraviolet and Visible Radiation model (TUVv5.2; https://www.acom.ucar.edu/Models/TUV/Interactive_TUV; last access: 27 January 2023) following the approach of recent studies (Lyu et al., 2019, 2016). Specifically, the geographical coordinates, date, and time were initialized into the TUV model to derive photolysis rates and solar radiation. We obtained the scaling factor by comparing the observed with modeled solar radiation and used this scaling factor to scale the TUV-model-derived photolysis rates. A dilution rate of $\mathrm{3}/\mathrm{86}\mathrm{400}$ s⁻¹ was applied for all non-constraint species and simulation days through a stepwise sensitivity test by adjusting it from $\mathrm{1}/\mathrm{86}\mathrm{400}$ to $\mathrm{5}/\mathrm{86}\mathrm{400}$ s⁻¹ (see details in Sect. S1) for the best reproduction of O₃. For each model run (i.e., each daily model simulation), it was performed on a daily basis with intervals of 24 h spanning from 00:00 to 23:00 LT (local time), and each individual model simulation was run to reach one-day diurnal steady state. The detailed descriptions of box model operation were provided in our previous study (Li et al., 2021).

Since the box model simulations are conducted with intervals of 24 h spanning from 00:00 to 23:00 LT (Wang et al., 2018), the entire campaign's observations were classified into four patterns of timescale (i.e., 5-month, monthly, weekly, and daily) as diurnal-average formats for model input (Fig. 1). Note that some days or weeks were not modeled due to some missing data in the measurements. Nevertheless, the total simulation number at the daily (i.e., 100, 81, and 114 d for TZ, BJ, and XD respectively) or weekly (i.e., 21, 20, and 19 weeks for TZ, BJ, and XD respectively) scale was representative of the 5-month campaign. Specifically, the entire campaign data classified as four patterns of timescale were modeled as base runs. Then we performed the sensitivity modeling to calculate the relative incremental reactivity (RIR) of precursors by adjusting the input concentrations in the base runs (see next section) (Lu et al., 2010a).

Figure 1 Schematic diagram of the dataset treatment to derive four patterns of timescale for 0-D box model input. Note that the four patterns (i.e., 5-month, monthly, weekly, and daily) were the diurnal average of the initial dataset. This diagram takes one site and several input measurements (temperature, toluene, and NO₂) as examples.

2.3 Calculation of net O_x production rate P(O_x) and relative incremental reactivity (RIR)

Considering the rapid chemical titration of NO to NO₂ in the presence of O₃, the concept of “total oxidant” (O_x= O₃+ NO₂) has been widely used to represent the actual photochemical production of O₃ (Lu et al., 2010a). Similar to those described in previous studies using the 0-D box model (He et al., 2019; Lyu et al., 2016), the net or in situ O_x production rate (P(O_x)) is defined as the difference between the O_x gross production rate (G(O_x)) and the O_x destruction rate (D(O_x)), which is formulated in accordance with Eq. (1):

$$\begin{array}{}\text{(1)}& P\left({\mathrm{O}}_x\right)=G\left({\mathrm{O}}_x\right)-D\left({\mathrm{O}}_x\right).\end{array}$$

The O_x gross production rate (G(O_x)), or the total chemical production of O_x, is calculated by summing the rates of oxidation of NO by HO₂ and RO₂ radicals in accordance with Eq. (2):

$$\begin{array}{}\text{(2)}& \begin{array}{rl}G\left({\mathrm{O}}_x\right)& =k_{{\mathrm{HO}}_{\mathrm{2}}+\mathrm{NO}}\left[{\mathrm{HO}}_{\mathrm{2}}\right]\left[\mathrm{NO}\right]\\ & +\sum k_{{\mathrm{RO}}_{\mathrm{2},i}+\mathrm{NO}}\left[{\mathrm{RO}}_{\mathrm{2},i}\right]\left[\mathrm{NO}\right],\end{array}\end{array}$$

The O_x destruction rate (D(O_x)), or total chemical loss of O_x, is calculated by summing O₃ photolysis, the reaction of O₃ with OH, HO₂, and alkenes, and the reaction between NO₂ and OH, as described by Eq. (3):

$$\begin{array}{}\text{(3)}& \begin{array}{rl}\mathrm{D}\left({\mathrm{O}}_x\right)& =k_{{\mathrm{O}}^{\mathrm{1}}\mathrm{D}+{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}\left[{\mathrm{O}}^{\mathrm{1}}\mathrm{D}\left]\right[{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\right]+k_{\mathrm{OH}+{\mathrm{O}}_{\mathrm{3}}}\left[\mathrm{OH}\right]\left[{\mathrm{O}}_{\mathrm{3}}\right]\\ & +k_{{\mathrm{HO}}_{\mathrm{2}}+{\mathrm{O}}_{\mathrm{3}}}\left[{\mathrm{HO}}_{\mathrm{2}}\right]\left[{\mathrm{O}}_{\mathrm{3}}\right]+k_{\mathrm{alkenes}+{\mathrm{O}}_{\mathrm{3}}}\left[\mathrm{alkenes}\right]\left[{\mathrm{O}}_{\mathrm{3}}\right]\\ & +k_{\mathrm{OH}+{\mathrm{NO}}_{\mathrm{2}}}\left[\mathrm{OH}\right]\left[{\mathrm{NO}}_{\mathrm{2}}\right].\end{array}\end{array}$$

Concentrations of radicals and intermediates are obtained from the outputs of the 0-D box model. The k values in Eqs. (2) and (3) represent the rate constants of the corresponding reactions, respectively. The subscript “i” in Eq. (2) represents the individual RO₂ species.

Additionally, relative incremental reactivity (RIR) has been widely used as a metric to quantify the O₃–precursor relationship, and it can be derived from the 0-D box model (MCMv3.3.1) by changing the input mixing ratios of its precursors (Sillman, 2010; Xue et al., 2014a). The RIR is defined as the ratio of percentage change in net O_x (O_x= O₃+ NO₂) production rate P(O_x) (Li et al., 2021) to the percentage change of the concentration of precursor X. The RIR of a specific precursor X is described in Eq. (4):

$$\begin{array}{}\text{(4)}& \mathrm{RIR}\left(X\right)={\displaystyle \frac{\left[P{\mathrm{O}}_{\mathrm{x}}\left(\mathrm{CX}\right)-P{\mathrm{O}}_x(\mathrm{CX}-\mathrm{\Delta }\mathrm{CX})\right]/P{\mathrm{O}}_{\mathrm{x}}\left(\mathrm{CX}\right)}{\mathrm{\Delta }\mathrm{CX}/\mathrm{CX}}}.\end{array}$$

Here, X is a specific precursor (i.e., NO_x, CO, or grouped or individual VOC species), CX is the measured concentration of precursor X, and ΔCX is the hypothetical concentration change (ΔCX $/$ CX = 10 % in this study in accordance with the previous studies (Lyu et al., 2016; Wang et al., 2018)). PO_x(CX) represents the simulated O_x production rate in a base run, whereas PO_x(CX–ΔCX) is the simulated O_x production in a second run with a hypothetical concentration change of species X. Obviously, a higher positive value of RIR(X) suggests a more effective way of reducing the ambient O₃ production rate by reducing X (Ling et al., 2011; Zhang et al., 2008a).

In this study, the O₃ precursors were divided into four major categories, including anthropogenic VOCs (AVOCs), biogenic VOCs (BVOCs, only isoprene in this study), CO, and NO_x (Tan et al., 2019). AVOCs were further divided into three subcategories: alkanes, aromatics, and alkenes^* (the asterisk denotes anthropogenic alkenes, excluding isoprene in this study; Yu et al., 2020a). As mentioned, RIR method was applied mainly to evaluate the O₃–NO_x–VOC sensitivity and to determine the photochemical regimes among four patterns of timescale. Thus, we calculated the RIR values of major precursor groups (i.e., AVOCs, BVOCs, CO, NO_x, alkanes, alkenes*, and aromatics) to further quantify the O₃–precursor relationship.

In general, O₃ formation chemistry is usually classified into three regimes (i.e., VOC-limited, transitional, and NO_x-limited; He et al., 2019; Wang et al., 2018). In this study, RIR${}_{{\mathrm{NO}}_x}$ $/$ RIR_AVOC (the ratio of two RIR values) was used as a metric to classify the photochemical regimes (Li et al., 2021). Specifically, a RIR${}_{{\mathrm{NO}}_x}$ $/$ RIR_AVOC value of less than 0.5 was defined as a VOC-limited regime; a value greater than 2 was defined as a NO_x-limited regime; and a value from 0.5 to 2 was defined as a transitional regime (see Sect. S2 and Table S4) (Li et al., 2021).

3 Results and discussion

3.1 Overview of the field campaign

Figure 2 shows the time series of measured meteorological parameters and O₃, as well as its precursors at the three sites during the whole campaign. In general, the temperature (T) and relative humidity (RH) were basically consistent at the three sites, while the wind speeds were different, which suggests that the three sites had an overall consistent meteorological condition. In addition, the time series of UV-A radiation was shown in Fig. 2d, which was only available from one urban site of Zibo but was expected to represent the whole Zibo city in this study. Following the protocol of the previous studies (Lyu et al., 2019; Y. Wang et al., 2017; Xue et al., 2014c), the time series of photolysis rates (e.g., $J_{{\mathrm{NO}}_{\mathrm{2}}}$ (Fig. 2e) and $J_{{\mathrm{O}}^{\mathrm{1}}\mathrm{D}}$ (Fig. 2f)) were calculated from TUVv5.2 model and further scaled from UV-A radiation measurement.

Figure 2 Time series of meteorological parameters, O₃, and its precursors (i.e., CO, NO_x, VOCs) throughout the whole campaign at the three sites in Zibo.

As shown in Fig. 2g, we found that severe O₃ pollution was observed at the three sites throughout the whole campaign. According to our measurements at the three sites in Zibo, the averaged O₃ concentration during the whole observational period was around 50 ppbv, while the daily maximum of O₃ concentrations for some extremely polluted periods were nearly 120–150 ppbv (Fig. 2g). Interestingly, the O₃ concentrations at the three sites were generally consistent, while the levels of its precursors (e.g., VOCs, NO_x) were obviously different (Fig. 2h–k), which implies the site-to-site variation of O₃ formation chemistry for the whole city of Zibo.

Generally, OH reactivity (or OH loss rate, k_OH) is widely applied to quantify the capacity of OH consumption by VOCs (Tan et al., 2019). According to Table S3, the BVOC reactivity (k_BVOC, 3.5±4.1 s⁻¹) in TZ was highest among the three sites. As BJ was mainly influenced by the emission from urban region, it showed the highest AVOC reactivity (k_AVOC, 6.8±6.3 s⁻¹) and NO_x level (31.1±28.6 ppbv). In addition, XD showed the highest level of alkenes* reactivity of 4.0±3.2 s⁻¹ within the three sites, and the local petrochemical industry nearby the XD area may explain such a characteristic (Li et al., 2021).

3.2 Evaluation of box model performance

The measured O₃ concentrations were not constrained in our MCMv3.3.1 box model calculation; thus the model performance could be quantitatively assessed by comparing the modeled O₃ (from base runs) with the measured O₃. Figures S3–S8 show the time series of simulated and observed O₃ concentrations at four patterns of timescale. In most cases, the box model simulation could accurately capture the level and variation trend of the observed O₃. However, on some days, the modeling results underestimated or overestimated the O₃ concentrations – in particular, nocturnal O₃ concentrations were underestimated. Such discrepancies between the simulated and observed O₃ were likely due to limitations in the explicit representations of atmospheric and transport processes (i.e., the horizontal and vertical transport processes of ground ozone) by the 0-D modeling approach (Lyu et al., 2019; Yu et al., 2020b). Specifically, ozone simulated by the 0-D box model is considered in the same way as in situ photochemical processes from its precursors. Unlike the 3-D air quality model, 0-D box modeling usually simplifies the representation of the physical processes (i.e., deposition and advection; Lu et al., 2010a; Sillman, 2010). Note that some adjustable parameters (e.g., radiation scheme, dilution rate) remained consistent in all of our model calculations, which ensured the comparability of model results to the greatest extent.

The index of agreement (IOA; Li et al., 2021; Lyu et al., 2016), Pearson's correlation coefficient (r), and root-mean-square error (RMSE) were jointly used as statistical metrics to quantify the goodness of fit between the simulated and observed O₃ concentrations. Table S5 summarizes these statistical metrics for each site at various patterns of timescale. Because any single statistical metric has its own limitations, using these three indicators conjointly provided a more comprehensive evaluation of the model performance (Su et al., 2018b). Generally, higher IOA and r and lower RMSE indicate better agreement between the simulated and observed values (Wang et al., 2018; Willmott, 1982). As shown in Table S5, slightly reduced correlation was observed as the timescale changed from the wider (i.e., 5-month scale) to the narrower (i.e., daily scale) pattern, which is understandable because of the enlarged statistical samples in the narrower pattern of timescale.

In summary, TZ showed the best performance in terms of the box model simulation, followed by XD and then BJ, regardless of any statistical metrics or different patterns of timescale, which may be associated with the optimized dilution rate for non-constraint species in model configuration. The overall model performance in this study (i.e., a day-to-day IOA of approximately 0.90 for TZ) was close to or slightly better than those reported in previous studies, such as IOA = 0.74 in Hong Kong (Liu et al., 2019), IOA = 0.74 in Wuhan (Lyu et al., 2016), and IOA = 0.90 in Jiangmen (He et al., 2019). According to the above evaluation of base runs, our modeled results were acceptable for the subsequent O₃–precursor relationship analysis described in the following sections.

3.3 Month to month

Figure 3a–b presents the monthly RIR values of the major precursor groups at each site, and the large variability of the O₃–precursor relationship at a spatiotemporal scale (i.e., site to site and month to month) was observed. Specifically, in most months, XD generally showed the highest RIR_AVOC among the three sites, followed by BJ and then TZ. In addition, RIR_BVOC showed a similar level to RIR_AVOC in TZ but much less than RIR_AVOC in BJ and XD, which can be explained by the observed higher BVOC reactivity in TZ than in the other two sites (see Fig. S9 and Table S3). Also, almost all the precursor groups showed positive RIR values, except for the negative RIR${}_{{\mathrm{NO}}_x}$ that appeared in BJ and XD in September. In addition, the RIR_CO values at the three sites suggested its limited role in O₃ formation at the three sites compared with other major categories of O₃ precursors. Among the three subcategories of AVOCs, alkenes* always had the highest RIR values, followed by aromatics, while the contribution of alkanes to O₃ formation can be ignored due to their near-zero RIR values. That sequence of O₃–AVOC sensitivity (alkenes^* > aromatics > alkanes) indicated by the RIR analysis was consistent with previous studies in some other Chinese cities (Su et al., 2018b; Tan et al., 2019). Significant monthly variations of O₃, NO_x, CO, VOC reactivity, and TVOC $/$ NO_x ratios (in ppbC ppbv⁻¹ as a widely used simple metric to determine the photochemical regime; National Research Council, 1991) were also observed from May to September (see Fig. S9 and Table S3) at the three sites. For example, the BVOC reactivity in TZ showed the highest level among the three sites during the whole campaign, and the AVOC reactivity in BJ showed more considerable variations in different months, which indicated spatial and temporal variations of local primary emissions for O₃ precursors in the city of Zibo.

Figure 3 Time series of month-to-month RIR values of major precursor groups and RIR${}_{{\mathrm{NO}}_x}$ $/$ RIR_AVOC at three sites (TZ, BJ, and XD) in Zibo. The dashed green line indicates RIR${}_{{\mathrm{NO}}_x}$ $/$ RIR_AVOC= 0.5 and 2.

Figure 3c shows monthly RIR${}_{{\mathrm{NO}}_x}$ $/$ RIR_AVOC at each site, which clearly reveals the spatial and temporal variations in photochemical regimes. For instance, the photochemical regime at the TZ site was considered to be a transitional regime in May, a NO_x-limited regime in June and July, and a VOC-limited regime in August and September; on the other hand, for a specific month like June, NO_x-limited, VOC-limited, and transitional regimes were generally identified for TZ, BJ, and XD, respectively. Figure 5b shows good consistency between monthly TVOC $/$ NO_x and RIR${}_{{\mathrm{NO}}_x}$ $/$ RIR_AVOC, suggesting that the changes of local emissions for O₃ precursors may partially explain the considerable variation of O₃ formation chemistry in different months.

Figure 4 Time series of week-to-week RIR values of major precursor groups and RIR${}_{{\mathrm{NO}}_x}$ $/$ RIR_AVOC at three sites (TZ, BJ, and XD) in Zibo. The blue lines in (g)–(i) are the three-points moving average of RIR${}_{{\mathrm{NO}}_x}$ $/$ RIR_AVOC values.

Figure 5 The correlations of TVOC $/$ NO_x with RIR${}_{{\mathrm{NO}}_x}$ $/$ RIR_AVOC at multiple patterns of timescale at the three sites in Zibo.

3.4 Week to week

Figure 4 shows the time series of week-to-week RIR values of major precursor groups and RIR${}_{{\mathrm{NO}}_x}$ $/$ RIR_AVOC at three sites in Zibo. Compared with month-to-month results, Fig. 4 further reveals the O₃–precursor relationship with more information on temporal trends. The temporal variations in weekly RIR_AVOC at the three sites generally decreased and then increased, whereas weekly RIR${}_{{\mathrm{NO}}_x}$ represented an opposite temporal variation during the entire campaign. Additionally, weekly RIR_BVOC showed a trend of first decreasing and then increasing at TZ, while it did not show clear temporal variations at BJ and XD due to low values (Fig. 4a–c). In general, RIR_alkanes, RIR${}_{{\mathrm{alkenes}}^{*}}$, and RIR_aromatics showed a tendency consistent with that of the RIR_AVOC at the three sites (Fig. 4d–f). Overall, these phenomena were consistent among the three sites, though the magnitude of RIR values varied site to site. In parallel, the temporal changing of O₃ precursors (e.g., AVOCs, NO_x) was also observed at the three sites during the entire campaign (see Fig. S10). For example, the weekly NO_x concentration showed an overall trend of first decreasing and then increasing, while the AVOC reactivity showed a different temporal variation. Given the moderate correlation between weekly TVOC $/$ NO_x and RIR_AVOC $/$ RIR${}_{{\mathrm{NO}}_x}$ (Fig. 5c), the temporal variations of RIR values and O₃ formation chemistry at the three sites may be partially elucidated by the emission changes of O₃ precursors.

As shown in Fig. 4g–i, all three sites showed similar temporal trends in terms of RIR${}_{{\mathrm{NO}}_x}$ $/$ RIR_AVOC, as it increased first and then decreased, though the magnitude of RIR${}_{{\mathrm{NO}}_x}$ $/$ RIR_AVOC varied largely at each site. Such site-to-site variability of RIR${}_{{\mathrm{NO}}_x}$ $/$ RIR_AVOC suggests that the photochemical regime at a local scale was mainly influenced by local emissions. By contrast, the site-to-site synchronization in the temporal trend of RIR${}_{{\mathrm{NO}}_x}$ $/$ RIR_AVOC suggests that the photochemical regime at a local scale may also be influenced by the emissions in a regional area. Therefore, the long-term, week-to-week RIR${}_{{\mathrm{NO}}_x}$ $/$ RIR_AVOC of multiple sites can further reflect the variability of the ozone formation regime at a large geographic scale.

3.5 Day to day

In this section, O₃–precursor relationship at the narrowest pattern of timescale was identified in detail. Figures S11–S12 show the time series of daily RIR values at three sites in Zibo, where the temporal trend of RIR values was consistent with that at a weekly scale (Fig. 4). Additionally, the time series of daily RIR${}_{{\mathrm{NO}}_x}$ $/$ RIR_AVOC (Fig. S13) showed more irregular variations in temporal trends during the entire campaign, though such temporal trends were overall consistent with that of the weekly scale in Fig. 4g–i. In summary, the time series of RIR values from the daily scale can provide more informative variations and characteristics of the O₃–precursor relationship in terms of temporal trends.

Table 2 summarizes the number of days and proportions that were classified into the three photochemical regimes across each site and each pattern of timescale. Near-consistent proportions of O₃ formation regimes (using RIR${}_{{\mathrm{NO}}_x}$ $/$ RIR_AVOC as a metric) were shown among multiple patterns of timescale, whereas a variability of proportion occurred among the three sites. The proportions of photochemical regimes changed according to the timescale that was varied from wider to narrower patterns. Taking TZ as an example, 20 % (monthly) and 26 % (daily) of the time was considered to be a VOC-limited regime. The number of days and proportions for the photochemical regimes summarized at four patterns of timescale can reveal a more plausible and comprehensive variation in ozone formation chemistry. Compared with patterns of monthly and weekly scales, the results derived at a daily scale can reveal the temporal variability of photochemical regimes more comprehensively. Note that the photochemical regime proportion obtained from the day-to-day scale has an advantage due to the large number of statistical samples.

Table 2 Summary of the number of days (for model calculation) and proportions that were classified into the three photochemical regimes across each site and multiple patterns of timescale.

3.6 Comparison among different patterns of timescale

This section gives a more comprehensive understanding of the campaign-averaging O₃–precursor relationship by comparing the similarities and differences of the results from various patterns of timescale. The overall O₃–precursor relationship for the entire campaign can be quantified by averaging the RIR values from the individual simulation runs depending on the chosen timescale (e.g., five simulation runs for monthly scale in this study). Therefore, four sets of logical and comparable results can be derived to represent the campaign-averaging O₃–precursor relationship, as four patterns of timescale (i.e., 5-month, monthly, weekly, and daily) were treated in this study.

Figure 6 Distribution of RIR values of major precursor groups in multiple patterns of timescale at three sites (TZ, BJ, and XD) in Zibo.

Figure 6 shows the averaged RIR values of the major precursor groups at different patterns of timescale. As the timescale changed from a wider (i.e., 5-month scale) to narrower (i.e., daily scale) pattern, all three sites showed increases in the means of RIR_AVOC and RIR${}_{{\mathrm{alkenes}}^{*}}$, as well as decreases in averaged RIR${}_{{\mathrm{NO}}_x}$, whereas the averaged RIR of other precursors (i.e., BVOCs, CO, alkanes, and aromatics) did not vary obviously (see Table S6). Comparing with the O₃–VOCs–NO_x sensitivity at the daily scale, the results obtained at the 5-month scale underestimated O₃–AVOCs sensitivity (indicated by averaged RIR values) by 48 % (TZ), 66 % (BJ), and 49 % (XD) and overestimated O₃–NO_x sensitivity by 37 % (TZ), 142 % (BJ), and 144 % (XD). We performed comprehensive uncertainty analysis for model input and output results, which was assessed through statistical methods (see details in Sect. 3.7). We found that the model-derived RIR values may become more uncertain when the input dataset was averaged into a wider diurnal pattern (i.e., 5-month scale), which may explain the discrepancy in RIR values between the 5-month scale and daily scale. We expect that such discrepancies derived from different patterns of timescale could widely exist in many other world areas. Note that the mean RIR values were generally consistent among the four patterns of timescale within a reasonable range (within the 25–75th quantile and standard deviation; see Fig. 6 and Table S4), suggesting that any selected pattern of timescale could reasonably derive the campaign-averaging O₃–precursor relationship.

Figure 7 Distribution of RIR${}_{{\mathrm{NO}}_x}$ $/$ RIR_AVOC (indicator of photochemical regime) in multiple patterns of timescale at three sites (TZ, BJ, and XD) in Zibo.

Figure 7 further shows the variations in photochemical regimes (defined by RIR${}_{{\mathrm{NO}}_x}$ $/$ RIR_AVOC; see Sect. S2 and Table S4 for details) for each pattern of timescale. Specifically, TZ was mainly considered to be a transitional regime for the entire campaign period, whereas its variations covered three photochemical regimes, which was consistent with the results from Table S6. BJ was generally identified as a VOC-limited regime, whereas some days were also grouped into a transitional regime. XD was considered to be primarily between a VOC-limited and transitional regime, and its variations also spanned three photochemical regimes. Compared with the 5-month pattern, it was further found that the averaged RIR${}_{{\mathrm{NO}}_x}$ $/$ RIR_AVOC from other timescale patterns (i.e., monthly, weekly, and daily) was higher (12 % to 20 % for TZ; 38 % to 153 % for XD) or lower (21 % to 65 % for BJ) than that from the 5-month scale. Note that the above discrepancies in photochemical regime derived from multiple patterns of timescale may influence the development of targeted O₃ control strategies. In summary, the photochemical regime derived by averaging RIR${}_{{\mathrm{NO}}_x}$ $/$ RIR_AVOC from the daily scale (see Table S6) suggests that the three sites mainly followed the sequence of TZ (1.34±1.39) > XD (0.67±1.49) > BJ (0.16±0.65).

In addition, the temporal variations of TVOC $/$ NO_x at different timescales were identified during the whole campaign, and good correlations between observed TVOC $/$ NO_x and model-derived RIR${}_{{\mathrm{NO}}_x}$ $/$ RIR_AVOC at four patterns of timescale were also found (see Fig. 5). Such consistency suggests that both metrics can reasonably reflect the variation of photochemical regimes, which can also improve the reliability of our box model simulation.

The consistency and difference of model output (summarized in Table S7) are quantified by the statistical methods of Pearson's correlation coefficient (Hu et al., 2018) and paired-samples t-test analysis (Wang et al., 2016). In particular, we assess and compare the degree of significance of the differences among multiple patterns of timescale by means of the p values (a statistical significance assuming that p < 0.05) through paired-samples t-tests and Wilcoxon matched-pairs signed-rank tests (non-parametric statistics; Chiclana et al., 2013). Figure 8a shows that high Pearson's correlation coefficients (with values all above 0.85, p < 0.01) were found among four patterns of timescale and that the higher correlation coefficient was identified between the two closer patterns. Figure 8b–c shows that the differences among multiple patterns of timescale were non-significant using paired-samples t-test analysis and Wilcoxon matched-pair signed-rank tests, respectively. Furthermore, their results indicate that a more significant difference was recognized between the two distant patterns (e.g., daily and 5-month), which is consistent with the results of Pearson's correlation analysis. It is noted that the discrepancy between the two distant patterns was not significant but non-negligible (e.g., p= 0.092 of the Wilcoxon matched-pairs signed-rank test between 5-month and daily patterns).

Figure 8 The statistical analysis results of RIR values (from Table S6) at multiple patterns of timescale: (a) Pearson's r correlation analysis (all the results have passed statistical significance, assumed to be p < 0.01); (b) paired-samples t-test analysis (^*p values refer to differences with a statistical significance assumed to be p < 0.05); (c) Wilcoxon matched-pairs signed-rank test (^*p values refer to differences with a statistical significance, assumed to be p < 0.05).

The influence of different patterns of timescale on deriving RIR values from individual AVOC species was further investigated. Briefly, quantifying the relative contribution of individual AVOCs to O₃ formation based on RIR calculation is beneficial to the development of cost-effective AVOC control strategies (Zhang et al., 2021). Figure 9 shows the averaged RIR values of individual AVOC species (i.e., top 10) at different patterns of timescale (i.e., 5-month, month to month, week to week) at three sites in Zibo. As shown in Fig. 9, the 10 individual AVOC species at the three sites were selected according to the top 10 highest RIR from the 5-month pattern. All three sites showed that the RIR of individual AVOC species increased gradually as the timescale changed from the wider (i.e., 5-month) to narrower (i.e., weekly) pattern, which was consistent with the earlier discussion (see Fig. 6 and Table S6) of O₃–AVOC sensitivity derived from four patterns of timescale. The results also indicate that the choice of timescale pattern has a limited effect on deriving high-ranking AVOC species (i.e., top 10) based on RIR calculations.

Figure 9 Averaged RIR values of individual AVOC species (top 10) at different patterns of timescale at three sites (TZ, BJ, and XD) in Zibo. The error bars represent the standard deviations of the mean.

3.7 Uncertainty analysis

The uncertainty of model input, which is embedded in a pre-processed dataset with multiple patterns of timescale, was quantified in this section. As shown in Fig. 1, the daily simulation used the individual daily pattern to constrain the model, while the input dataset of averaged diurnal patterns (i.e., weekly, monthly, and 5-month) is treated by averaging the individual daily pattern into different timescales. This averaging approach will conceal the temporal variations of O₃ precursors and meteorological factors, particularly for a long-term observational campaign. Figure S14 shows the distributions of the standard deviations for OH reactivity (k_OH) or the concentration of O₃ precursor groups at three averaged patterns of timescale at the three sites. As the timescale changed from a wider (i.e., 5-month scale) to narrower (i.e., weekly scale) pattern, the uncertainty (indicated by the average, median, and 25 %–75 % quantile) decreased accordingly. In addition, meteorological factors such as temperature and irradiation also play an important role in O₃ formation, and it is specifically noted that these meteorological parameters can vary greatly over a long observational period (Boleti et al., 2020; Liu et al., 2019; Weng et al., 2022). Therefore, the masked temporal variation of these meteorological factors behind the averaged input dataset would also result in model uncertainty.

Moreover, it has been widely recognized that the uncertainty for 0-D box model simulation mainly arises from the constraint of observation datasets and the configuration of model schemes. Note that constraints with more species from measurements (or including as many species as possible) would lower its uncertainty from the chemical box model simulation (Wolfe et al., 2011, 2016). Nevertheless, due to the measurement limitation in our field campaign, we are unable to measure some important atmospheric species (i.e., HONO and oxygenated VOCs (OVOCs)), and these may raise uncertainty in the box model simulation. For instance, Xue et al. (2021) performed a sensitivity test for HONO constraint in their box model simulation, and they showed that no HONO constraint would lead to the O₃ photochemical production rate decreasing by 42 %. More recently, Wang et al. (2022) obtained a comprehensive VOC dataset at Guangzhou, and their results showed that box model simulation without OVOC constraints would underestimate the productions of RO_x and O₃. Besides, both gaseous HNO₃ and organic nitrates can result in interferences to NO_x measurements by means of the chemiluminescence technique, which may raise uncertainty in our box modeling (Ge et al., 2022; Uno et al., 2017; Xu et al., 2013). Since the accurate NO_x measurement is essential in determining the photochemical regime, more in-depth studies on NO_x measurement uncertainty in box model simulations are required in the future. In addition, the parameter configuration of model schemes is essential to derive a reliable and valid model output, such as dilution rate as an important technical model parameter. We performed a stepwise sensitivity test for this parameter to obtain an optimized dilution rate and assigned it to all non-constraint species, which can reduce uncertainty in box model simulation (see details in Sect. S1). Also, the dry and/or wet deposition of pollutants is an important atmospheric physical process, which has been mostly parameterized in emission-based chemical transport modeling but is very limited in box modeling, as most of the primarily emitted species are already constrained from measurements. Xue et al. (2014c) considered O₃ deposition in box model simulations, and their results showed a negligible contribution of O₃ deposition to total O₃ destruction rates. As for this work, we are unable to consider the deposition due to the difficulty in representing and parameterizing this term in the 0-D box model. Nevertheless, deposition of O₃ and other species may be one of the uncertainties during box model simulation, which is worth further study in the future.

4 Summary and implications

Our present results suggest that comprehensively understanding multiple patterns of timescale is conductive to formulating a more accurate and robust O₃ control strategy. Specifically, as identified from the narrower patterns of timescale (i.e., weekly and daily), the site-to-site photochemical regime indicated by RIR${}_{{\mathrm{NO}}_x}$ $/$ RIR_AVOC showed various magnitudes but a synchronous temporal trend. This indicates that the O₃ formation regime in a city area can be influenced by local and regional emissions jointly. The reason behind this phenomenon is not clear at present, and we believe that further investigation on the synergetic effect of local and regional emission reduction for O₃ control would help elucidate this observation. It was also found that the campaign-averaging photochemical regimes showed overall consistency but non-negligible variability among the four patterns of timescale, which was mainly due to the embedded uncertainty in the model input dataset with averaged diurnal patterns. This implies that comparison among multiple patterns of timescale based on RIR analysis is useful to derive the O₃–precursor relationship more accurately and reliably.

Moreover, the high-ranking AVOC species (i.e., top 10) based on RIR calculations were overall consistent between the narrow to wide patterns of timescale. Table S8 summarizes the total run number of box models for different patterns of timescale. It is known that large-scale computing capacity and computational efficiency were required in the narrower pattern of timescale (e.g., 2760 simulation runs at the weekly scale in this study). Considering the difficulties of performing long-term and continuous online measurements in some environments, it is also advisable to identify the high-ranking VOC species from the campaign-averaging diurnal pattern in box model simulation.

In this study, we explored the non-linearity of the O₃–precursor relationship in a way driven by the actual daily, weekly, and monthly variability around the distribution. Our results highlight the importance of quantitatively testing the impact of different timescales on photochemical regime determination, as there is uncertainty embedded in the model input dataset when averaging individual daily pattern into different timescales. Such understanding would be complementary in developing more accurate O₃ pollution control strategies, particularly as the long-term O₃–precursor observations (e.g., from several months to years) are becoming more available than before in many places throughout China. In addition, site-to-site differences of model-derived photochemical regimes also underline the importance of developing targeted O₃ control strategies for different areas on a city scale. Specifically, according to the averaged RIR${}_{{\mathrm{NO}}_x}$ $/$ RIR_AVOC of the daily pattern, the derived photochemical regime was transitional for TZ (suburban) and XD (suburban), while it was VOC-limited for BJ (urban). This implies that, for mitigating ozone pollution in the city of Zibo, more endeavors should be devoted to the anthropogenic VOC reduction in urban areas while strengthening the synergetic mitigation of VOC and NO_x emissions at the same time in other suburban areas. Although the above implications for O₃ control were derived from a case study in a major prefecture-level city (Zibo) of northern China, the approach of integrating multiple patterns of timescale developed in the present work can be used in other regions, particularly in relation to the ongoing “One City One Policy” campaign (2021–2023) for O₃ control in many cities in China.