Observed changes in the temperature dependence response of surface ozone under NOx reductions

Due to the strong temperature dependence of surface ozone concentrations (O3), future warmer conditions may worsen ozone pollution levels despite continued efforts on emission controls of ozone precursors. Using long-term measurements of hourly O3 concentrations co-located with NOx concentrations in stations distributed throughout Germany, we assess changes in the climate penalty, defined as the slope of ozone-temperature relationship during the period 1999-2018. We find a stronger temperature sensitivity in the urban stations over the southwestern regions, especially in the first period of the study 5 (1999-2008). We show a decrease in the climate penalty in most of stations during the second period of the study (2009-2018), with some exceptions (e.g. Berlin) where the climate penalty did not show significant changes. A key motivation of this study is to provide further insights into the impacts of NOx reductions in the O3-temperature relationship. For that, we propose a statistical approach based on generalized additive models (GAMs) to describe ozone production rates, inferred from hourly observations, as a function of NOx and temperature, among other variables relevant during the O3 production. We find lower 10 O3 production rates during the second period (2009-2018) at most stations and a decreasing sensitivity to temperature, pointing out that lowering NOx concentrations resulted in decreasing O3 production rates. However, we also observe changes in the shape of the function representing the O3-temperature relationship, which indicate that NOx reductions alone can not explain the changes in the temperature dependence of O3. Our analysis would suggest that decreasing NOx concentrations are not the only factor causing the observed changes in the climate penalty factor. 15

of the monitoring stations with co-located data (O 3 and NO x ) was based on the station type (background), station type area (rural, urban, suburban) and altitude (<1000m). Only the stations reporting more than 75 % of valid data out of all the possible data in each summertime were included in the study. We use the stations with at least 19 years with hourly co-located data within the whole period of study defined from 1999-2018. Here, summertime is referred to July-August-September (JAS), 5 with a strong O 3 -temperature relationship, particularly in Central Europe (Otero et al., 2018). A total of 29 stations meet the pre-processing criteria: 15 rural, 12 urban and 2 suburban stations. Despite that the spatial distribution of the measurement The meteorology was extracted from the ERA5 (Herbach and Dee, 2016), the latest climate reanalysis produced by the European Centre for Medium Range Weather Forecast (ECMWF) that provides hourly data on regular latitude-longitude 0.25ºx0.25º spatial resolution. The variables included in the analysis are air surface 2m-temperature ( • C), 10m u and v-component of wind (m s −1 ), boundary layer height (m) and relative humidity at 1000 hPa (%).

Climate penalty factor
A number of definitions have been used in the literature to characterise the ozone climate penalty, usually represented as the linear relationship between O 3 and temperature. Climate penalty values are normally computed using daily maximum summertime O 3 observations (1 or 8h average) and daily maximum temperature, although there is no standard definition 20 (Pusede et al., 2015). Here, we adopted one of the most common metric to represent the climate penalty (hereinafter, m O3T ) as the slope of the best fit line between long-term MDA8 concentrations and daily maximum temperature (Bloomer et al., 2009;Steiner et al., 2010;Otero et al., 2018). We first calculated the m O3T with a linear regression model applied separately for each station and each period (1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018). The general equation for the linear model can be written as follows: (1) 25 Then, we examined the difference between the slopes obtained for each period, introducing an interaction factor term in the linear model to quantify the slope differences: Y(t), T(t) are the time series of MDA8 and daily maximum temperature (respectively) for the whole period 1999-2018 and P is categorical variable with two categories: one representing the period 1999-2008 and another representing 2009-2018. 30 4 https://doi.org/10.5194/acp-2020-691 Preprint. Discussion started: 18 August 2020 c Author(s) 2020. CC BY 4.0 License.

Approximation of O 3 production rates from observations
Most of the previous works have used numerical models (Steiner et al., 2006), box model (Coates et al. 2016, plume model (LaFranchi et al., 2011) or analytical models (Pusede et al., 2014;Romer et al., 2018) to analyse the temperature-dependent mechanisms affecting the O 3 production. Here, we propose a new approach based on GAMs to examine changes in the O 3 production. We approximate the later by the rate of change of hourly O 3 concentrations as: The general O 3 budget equation can be expressed as: P O 3chem represents the chemical O 3 production rate, LO 3chem is the chemical loss rate and the last term MD represents the dynamical processes that influence O 3 concentrations, including mixing and dry deposition processes. These individual processes 10 can vary in strength and by location throughout the day.
As we aim to assess how NO x reductions influence the sensitivity of O 3 to temperature, we restrict our analysis to a time interval with an intense photochemical activity, which usually coincides with higher O 3 concentrations and warmer temperatures. Thus, the data was filtered to avoid including non-related photochemical processes that might mask the photochemistry in the daily O 3 production. First, we selected data after sunrise and until O 3 reaches the daily maximum value (usually in 15 the afternoon). In order to exclude some maximum values that might occur late in the afternoon or evening and are mostly related to prevailing meteorological conditions and transport processes (Kulkarni et al., 1993), the time was restricted to 17:00 H (local time). Then, a wind speed condition was used to exclude the hourly data when wind speed was higher than 3.2 ms −1 , which is the threshold value usually applied to define stagnant conditions (Horton et al., 2014). After a first inspection of the data, we found considerable differences in the minimum of NO x concentrations across some stations and periods, likely due 20 to the detection methods. To better establish a comparison between stations and periods, we applied an additional filter to remove NO x values below 5 µ gm −3 . The number of observations that met these conditions varies with each station type and on average a 20% (urban), 14%(rural), 18%(suburban) of the total data was used.

Modeling O 3 production rates with GAMs
GAMs (Hastie and Tibshinari, 1990;Wood, 2006) were used to examine variations in ∆O 3 (t) over the last two decades and 25 the changes in the relationship NO x -temperature given the observed downward trends of the O 3 sensitivity to temperature in the two periods of study 1999-2008 and 2009-2018. GAMs are useful tools for estimating non-parametric relationships whilst retaining clarity of interpretation (Wood, 2006). The relationship between the explanatory variables (henceforth covariates) and the response is described by smooth curves (splines, or potentially other smoothers). Such models have proven useful for studying the complex non-linear relationships between atmospheric chemical species and meteorological parameters (Carslaw 30 et al., 2007;Jackson et al., 2009;Barmpadimos et al., 2011;Boleti et al., 2019). GAMs allow including nonlinear interactions between covariates with different smoothers assumed for each covariate. The GAM can be formally written as: where g is the link function, X n are the explanatory variables and f n are the non-parametric smoothing functions. β 0 is the intercept and is an error term. If the response can be assumed to be normally distributed, the canonical link function is the identity. After a closer inspection of the residuals at the individual sites, we found non normally distributed residuals with problems in the tails. Alternatively, we used a scaled t distribution recommended for heavy tailed response variables (Wood et al., 2016). Thin plate regression was used as smoothers to describe a nonlinear relationship between the response and 2 5 covariates (interaction) (Wood, 2006). The smoothness of each function is controlled by the number of knots or effective number of degrees of freedom. Here, the smoothing parameters were estimated by restricted maximum likelihood (REML) (Wood, 2006).
The challenge in building a model that captures a lage proportion of the variability of ∆O 3 is to select the key covariates out of a large number of potential variables. As stated in the previous section, changes in O 3 concentrations depend on local pro-10 duction, involving many chemical reactions that vary with temperature, loss mechanisms that are sensitive to meteorological conditions and transport processes. Therefore, we chose the variables that are expected to have a major influence on O 3 production (e.g. NO x ). The photochemical nature of O 3 production is strongly influenced by temperature. In particular, temperature Daytime variation in the boundary layer height (BLH) significantly contributes to changes in O 3 production rates that tend to increase with a deepening BLH during sunny and warm days. (Haman et al., 2014). In addition to chemical and mixing processes, changes in O 3 concentrations are influenced by deposition. Therefore, additional covariates are the percentage of change of the boundary layer height growth rate (∆BLH) (in %) accounting for mixing processes, and vapour pressure deficit 20 (VPD) as it has been recognised as a key variable for dry deposition (Kavassalis and Murphy, 2017;Otero et al., 2018). The VPD was calculated from the corresponding hourly data of air temperature and relative humidity. Moreover, we included the O 3 concentrations from the previous hour (C O3 (t − 1)) and the MDA8 concentrations from the previous day (C M DA8 (t−24)) to represent the persistence of previous chemical conditions, (Pusede et al., 2015). We first started with a baseline model that included the nonlinear relationship between NO x and temperature as follows: f (T, NO x ) represents the interaction between temperature (T) and NO x concentrations and it is included as a tensor product (Wood, 2017). Observing the skewness of the NO x data led us to introduce a modification in the baseline model using a log transformation of NO x . Since we aim to build a parsimonious model to better explain the variability of ∆O 3 , we gradually added the covariates to the baseline model through a selection procedure. During the stepwise process, we also allowed inter-

10
For that, we examine time series of the annual 5th, 50th, and 95th percentiles calculated from daily NO x concentrations, assessing the trends (Kendall, 1975) and estimating its slope (Theil, 1950;Sen, 1968). Figure 2 shows annual 5th, 50th, and 95th percentiles calculated from daily NO x concentrations at some example stations located in Berlin, Rhineland-Palatinate and Saxony that are representative for each station type area and will be used below to present the modelling results. The They found larger trends in m O3T at high and moderate polluted clusters ant they argued that it might be due to NO x reductions.
Here, we found a general decrease in m O3T obtained from long-term data across different environments (i.e. rural, urban and suburban). Our results also pointed out significant differences in the m O3T across stations, with some polluted areas where the m O3T did not show significant changes with time (e.g. Berlin). As stated in the introduction, mechanisms controlling m O3T are not well established. A priori it is not evident what the impact of NO x reductions is in the O 3 sensitivity to temperature, in particular in rural environments. Therefore, we next examine the variability of ∆O 3 as a function of temperature and NO x in order to provide further insights into the nonlinear temperature-dependence of NO x and the potential impacts on the observed 5 m O3T .

Model performance
A final model including three interaction terms was designed from the selection procedure as the best fit to capture the ∆O 3 variability at most stations and periods. The model selection indicated that as variables were added and the model complexity increased (i.e. more interactions), the AIC decreased and the deviance explained increased (Fig. S4). The performance of the 10 GAMs was assessed by the adjusted r-squared for the model (R 2 ), defined as the proportion of the variance explained (Fig. 4).
The results showed similar R 2 values in both periods over most of the stations, with some exceptions where GAM-P1 seem to perform better than GAMP-P2 (e.g. over the region of Hessen). In general, GAMs showed a better performance over urban and suburban stations and 40% of the ∆O 3 variability was captured. The models performed poorly when applied to rural stations, they showed lower values of R 2 . This likely reflects that GAMs designed with the underlying assumptions of the interactions 15 between the selected covariates is better suited for urban and suburban areas than for rural regions.

Model interactions
Our approach is built upon a conceptual model (4)   Berlin NO x declined only by 7.5%.
We examine the ∆O 3 response to temperature under the mean NO x conditions for each period using GAM-P1 and GAM-P2 along with the prediction obtained from GAM-P1 that projects the ∆O 3 response in the second period 2009-2018 (prediction line in Fig. 5, rigth). In Berlin, the ∆O 3 response to temperature shows a similar increase with temperature in both periods.

15
In this case, the GAM-P1 prediction for the second period 2009-2018 is in a good agreement with the shape obtained from GAM-P2, which suggest that a decreasing temperature sensitivity of ∆O 3 could be explained by NO x reductions. The increase of ∆O 3 with temperature is also depicted in Rhineland-Palatinate. The prediction from GAM-P1 for the second period 2009-2018 reveals discrepancies at higher temperatures when comparing to the ∆O 3 response from GAM-P2. It can be noted that the prediction from GAM-P1 for the second period (prediction line, Fig.5) does not capture the steepness at temperatures above 20 20 • C showed by GAM-P2. Contrasting to the results in Berlin, the changes in the shape that represents the ∆O 3 as a function of temperature suggest that the NO x reductions would only partially explain the observed changes in the O 3 -temperature relationship, but rather an underlying effect is likely to influence the ∆O 3 at higher temperatures. We found similar features in the rest of the urban stations than in the example stations, with consistent interaction surfaces in terms of the ∆O 3 response to NO x and the temperature dependence (Fig. S5). As in Rhineland-Palatinate, the regression lines were slightly different when 25 comparing GAM-P2 and the projected ∆O 3 response under NO x reductions (Fig. S6), which reinforce our hypothesis of an underlying factor influencing the ∆O 3 -temperature relationship.
We further assess the effect of the temperature and NO x on ∆O 3 separately with GAM-P1 and GAM-P2 under fixed NO x and temperature conditions determined as the 10th, 50th and 90th percentiles of the corresponding distributions over the whole period of study (1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018). In contrast to the contour plots (Fig. 5), we now include the intercept and a constant value 30 (i.e. median) for the rest of the covariates, in order to further examine the summed effects. Table 1 summarizes the values of the covariates for the selected stations. Figure 6 shows ∆O 3 as a function of temperature. ∆O 3 estimates are generally  of NO x . We see similarities between the rural and urban stations in Berlin, in terms of the shape of the nonlinear relationship between temperature and NO x , which is expected due to the proximity between both stations ( Fig. 1). In contrast to the urban stations, the contours obtained from GAM-P1 and GAM-P2 are significantly different, particularly in Rhineland-Palatinate.

10
The ∆O 3 as a function of temperature under NO x mean conditions is also shown in Fig.8 (rigth). In Berlin, NO x concentrations declined by 28.8%, which could partially explain the decrease in ∆O 3 estimates during the second period 2009-2018.
However the shapes of the regression lines obtained from the GAM-P2 and the projected ∆O 3 response from GAM-P1 differ at temperatures below 20 • C (blue and green lines, Fig. 8). In Rhineland-Palatinate the temperature dependence is considerably lower than in Berlin and a flat regression line is shown by GAM-P2 for the second period with a 37% decrease of NO x 15 concentrations. The discrepancies between the projected ∆O 3 response and the ∆O 3 estimates from GAM-P2 are higher at temperatures below~20 • C. Overall, we found a larger variability among the rest of the rural stations considered in the study, in terms of the interaction surfaces NO x -temperature (Fig. S7). This is also reflected in the estimated temperature response of ∆O 3 when comparing GAM-P2 and the projected response using GAM-P1 (Fig.S8). increases of NO x > 7 µ gm −3 indicate a major decreased of ∆O 3 at medium (50th) and high (90th) temperatures in the second period.
Only two suburban stations were included in this study, in Berlin (DEBE051) and in Saxony (DESN045), both eastward located. The contours obtained in each period and station showed similar patterns than those found for urban stations, specially in Berlin (Fig. S9). The GAMs consistently reproduce the temperature dependence of ∆O 3 at higher temperatures and the 30 differences between the GAM-P2 and the projected ∆O 3 response to temperature with GAM-P1 were more evident in Saxony (Fig. S9, rigth). Figure 11 reveals the nonlinear relationship between VPD and the C O3 (t − 1) at the selected urban stations in Berlin and in Rhineland-Palatinate. In general, ∆O 3 tends to increase with higher levels of VPD (i.e. drier conditions) and low O 3 concentrations from the previous hour in both locations and periods. In the first period, the contribution of the interaction between VPD and persistent O 3 concentrations is similar at both locations, and the model shows maximum ∆O 3 at C O3 (t − 1) < 30 µ gm −3 and VPD > 0.70 kPa. In Berlin, the results obtained from GAM-P2 suggest that higher levels of VPD and low C O3 (t − 1) (~30 µ gm −3 ) lead to an increase of ∆O 3 , but the ∆O 3 tends to decrease faster with high C O3 (t − 1) concentrations (above is consistent with the patterns found in the urban and rural stations and ∆O 3 increases with higher VPD and low C O3 (t − 1) 15 concentrations (Fig. S10).
VPD is crucial and controls the stomatal conductance. Its effects can be summarised as follows: under high VPD levels (associated with high temperatures), plants cannot extract sufficient moisture from dry soils to satisfy the atmospheric demand for evapotranspiration (Teuling, 2018). In this situation of drought stress, plants close their stomata to reduce water loss and limit the uptake of ozone by vegetation. 20 Our results illustrate that the combination of high VPD and lower C O3 (t − 1) concentrations result in higher ∆O 3 (thus, less uptake of O 3 ). Moreover, given that O 3 concentrations are typically lower in urban environments due to the local scavenge of O 3 (NO titration), a larger contribution of the interaction of VPD and C O3 (t − 1) to ∆O 3 in the urban and suburban stations than in the rural stations is expected.

∆BLH and MDA8 from the previous day (C
The effect of mixing processes was introduced in the GAMs through the ∆BLH and C M DA8 (t−24). Figures 13 and 14 depict the interaction surfaces between the covariates ∆BLH and C M DA8 (t − 24) at the selected urban and rural stations in Berlin and Rhineland-Palatinate. In general, ∆O 3 is mainly dependent on changes in ∆BLH and it increases at higher ∆BLH, while the influence of C M DA8 (t − 24) on ∆O 3 is negligible for ∆BLH~< 30%. The results obtained from most of the stations at different environments (i.e. urban and rural) showed consistent shapes with the patterns described for the selected stations (not These interaction surfaces can be used to interpret the nonlinear relationship between ∆BLH and C M DA8 (t − 24) concentrations. As BLH grows, air is entrained from layers aloft and O 3 production rates can increase or decrase depending on the O 3 concentrations in this residual layer (Haman et al., 2014). We show that a rapid development of the BLH along with high C M DA8 (t−24) (from the previous day), likely stored at the residual layer, lead to an increase of ∆O 3 . Note that C M DA8 (t−24) concentrations seems to have an influence on ∆O 3 when the BLH rapidly changes. The effect of this interaction was slightly larger in most of the urban and suburban stations as compared to the rural stations, while small differences are observed when 5 comparing the patterns obtained from each period.

Summary and conclusions
We have examined the long-term O 3 sensitivity to temperature, as well as the modulation of this sensitivity by NO x concentrations, in a total of 29 stations over Germany during the period 1999-2018. Consistent with previous work, O 3 tends to increase strongly with temperature under high NO x conditions due to increased in-situ photochemical production, while lower levels of 10 NO x leads to a reduced O 3 sensitivity to temperature. Also consistent with previous work, we see a decreasing sensitivity of O 3 to temperature over our study period, coinciding with decreasing trends in NOx concentration.
In order to explain the trends in photochemical ozone production over our study period, we divided this period into two halves (1999-2008 and 2009-2018) and constructed sets of Generalized Additive Models (GAMs) based on hourly station observations of ozone and NO x concentrations, along with temperature, vapor pressure deficit, and boundary layer height from 15 a reanalysis product. We modeled the daily increase in O 3 concentrations as a function of three interaction terms accounting for phochemical production (dependent on NO x and temperature), dry deposition (dependent on vapor pressure defict and ozone concentrations from the previous hour) and mixing processes (dependent on the boundary layer height growth rate, and ozone concentrations from the previous day).
In most of the stations, the effect of the interaction term NO x -temperature was larger in the first period than in the second We conclude that NO x reductions have had an influence in the decreasing temperature sensitivity of O 3 , as shown in the GAMs for the second period 2009-2018, but that such reductions alone can not explain the changes in the observed O 3temperature relationship. We interpret these discrepancies as an underlying effect influencing the ∆O 3 that has not been included in the model.

35
The temperature-dependence of bioginec VOC emissions is well-known. In particular, biogenic isoprene emissions have a strong temperature dependence with critical implications on O 3 production, mostly during warmer summer days. Therefore, one plausible explanation for the changes in the shapes of the ∆O 3 -temperature relationship might be attributed to the accompanying effect of changes in biogenic emissions (along with NO x ) that are likely influencing the temperature-dependence of ∆O 3 and consequently the m O3T . We have shown a general decrease of ∆O 3 at higher temperatures, which may suggest that 5 enhanced temperature-driven biogenic emissions can result in ∆O 3 being more dependent of NO x (NO x -limited) and thus, more sensitive to the NO x controls. Our results have important implications for the implementation of mitigation strategies, specially when considering the effects of a warming climate. We expect that the methodology described herein can be applied to other locations with available long-term measurements to assess how NO x reductions have influenced the temperature dependence of O 3 . Consistent with previous work, we may anticipate that our approach will show changes in the climate penalty factor as well as in the sensitivity of ∆O 3 with temperature. Further analysis to examine in more detail the effect of NO x reductions in particular locations should be required.
In summary, the sensitivity of O 3 to temperature has decreased during the last period (2009-2018) over a great number 5 of the German stations considered in the study, including rural areas. Even though NO x reductions accomplished during the last decades have partially counteracted the O 3 climate penalty, our study highlights the relevance of considering the influence of additional factors controlling the O 3 -temperature relationship. Since observations of long-term dataset of VOCs are lacking, further analysis including short-term measurements of a suite of VOCs would be definitively required to quantify their contribution to the observed changes in the climate penalty.