Climate change penalty and benefit on surface ozone: a global perspective based on CMIP6 earth system models

This work presents an analysis of the effect of climate change on surface ozone discussing the related penalties and benefits around the globe from the global modelling perspective based on simulations with five CMIP6 (Coupled Model Intercomparison Project Phase 6) Earth System Models. As part of AerChemMIP (Aerosol Chemistry Model Intercomparison Project) all models conducted simulation experiments considering future climate (ssp370SST) and present-day climate (ssp370pdSST) under the same future emissions trajectory (SSP3-7.0). A multi-model global average climate change benefit on surface ozone of −0.96 ± 0.07 ppbv °C−1 is calculated which is mainly linked to the dominating role of enhanced ozone destruction with higher water vapour abundances under a warmer climate. Over regions remote from pollution sources, there is a robust decline in mean surface ozone concentration on an annual basis as well as for boreal winter and summer varying spatially from −0.2 to −2 ppbv °C−1, with strongest decline over tropical oceanic regions. The implication is that over regions remote from pollution sources (except over the Arctic) there is a consistent climate change benefit for baseline ozone due to global warming. However, ozone increases over regions close to anthropogenic pollution sources or close to enhanced natural biogenic volatile organic compounds emission sources with a rate ranging regionally from 0.2 to 2 ppbv C−1, implying a regional surface ozone penalty due to global warming. Overall, the future climate change enhances the efficiency of precursor emissions to generate surface ozone in polluted regions and thus the magnitude of this effect depends on the regional emission changes considered in this study within the SSP3_7.0 scenario. The comparison of the climate change impact effect on surface ozone versus the combined effect of climate and emission changes indicates the dominant role of precursor emission changes in projecting surface ozone concentrations under future climate change scenarios.


Introduction
Surface ozone is widely recognized as a pollutant having an impact on human health, crops, and ecosystems while tropospheric ozone plays a key role in the oxidizing capacity of the troposphere and is an important greenhouse gas (Monks et al 2015, Young et al 2018. Given the multifaceted dependence of tropospheric ozone on several meteorological parameters and various gaseous precursors, and its varying lifetime (from a few hours in polluted urban regions up to several months in the upper troposphere), an estimation of the climate change impact on future tropospheric and surface ozone levels is a complicated task with a regional context.
The future tropospheric and surface ozone changes at the global scale depend on changes in chemical ozone production and loss processes, stratosphere-troposphere exchange (STE), and deposition (Young et al 2013). The future changes in chemical ozone production and loss are strongly related to changes in anthropogenic and natural emissions of ozone precursors, but also to changes of meteorological parameters, such as temperature, water vapour and radiation. During the 21st century, changes in climate, ozone depleting substances (ODSs), and emissions of ozone precursor species are expected to be the major factors governing ozone concentration distribution in the stratosphere, the free troposphere, and at the surface (Fiore et al 2015, Revell et al 2015.
Focusing solely on the response of surface ozone to climate-driven changes (instead of the response to anthropogenic emission-driven changes of ozone precursor species), the 5th IPCC Assessment Report (Kirtman et al 2013) assessed that there is high confidence that in unpolluted regions, higher water vapour abundances and temperatures in a warmer climate enhance ozone destruction leading to lower baseline ozone levels, while there is medium confidence that in polluted regions surface ozone is expected to increase. Baseline ozone is defined as the observed ozone at a site when it is not influenced by recent, locally emitted or anthropogenically produced pollution (Jaffe et al 2018). The above-mentioned IPCC AR5 assessment is also supported by several recent studies based on model projections either with global or regional models. The majority of the regional studies are focusing on North America (Gonzalez-Abraham et al 2015, Val Martin et al 2015, Schnell et al 2016, He et al 2018, Nolte et al 2018, Rieder et al 2018, Gao et al 2020 or Europe (Colette et al 2015, Lacressonnière et al 2016, Schnell et al 2016, Fortems-Cheiney et al 2017, but there are also a few studies focused on East Asia (Lee et al 2015, Schnell et al 2016 and India (Pommier et al 2018). Over polluted regions of the world, such as in North America, Europe, and East Asia, model studies project a general increase of surface ozone levels in a future warmer climate, particularly during summertime as has been recently reviewed (Fu and Tian 2019). This future surface ozone increase, due to a warming climate in the absence of changes in anthropogenic polluting activities, reflects a climate change penalty on ozone (Wu et al 2008, Fiore et al 2015, Fu and Tian 2019.
However, the response of surface ozone to climate induced Earth system changes is complex due to counteracting effects of various processes affected by a warmer climate. Such processes, including STE, biosphere interactions, lighting and soil NO x emissions, wildfires, as well as anticyclonic stagnation conditions and heatwaves, can affect and modify future baseline and regional/local surface ozone levels (Fiore et al 2012, 2015, Fu and Tian 2019. Studies considering the individual effects of climate driven changes in specific precursor emissions or processes show increases in surface ozone under warmer climate for certain processes. This is indeed the case for enhanced STE and stratospheric ozone recovery (Kawase et al 2011, Sekiya and Sudo 2014, Hess et al 2015, Banerjee et al 2016, Meul et al 2018, Morgenstern et al 2018, Akritidis et al 2019, and the increase of soil NO x emissions (Wu et al 2008, Romer et al 2018, which can each lead to 1 and 2 ppbv increase in surface ozone.
Other temperature-dependent processes, including biosphere interactions affecting biogenic volatile organic compounds (BVOCs), natural methane emis-  (Lin et al 2017) and anticyclonic stagnation conditions (Jing et al 2017, Schnell and Prather 2017, Garrido-Perez et al 2018 are expected to play a key role in future surface ozone and even occurrence of high pollution events (e.g. in the case of wildfires and anticyclonic stagnation conditions) but their effects are difficult to quantify in isolation. For example, models can disagree on both the sign and magnitude of changes in surface ozone due to the individual model's climate sensitivity and varying level of complexity in the implementation of processes like temperature and land-use/land-cover sensitive BVOC emissions, deposition, CO 2 inhibition of isoprene emissions and model assumptions on the yields of isoprene nitrates (Squire et al 2015, Val Martin et al 2015, Schnell et al 2016, Pommier et al 2018. The direction of change of lightning activity in the future remains also highly uncertain (Clark et al 2017, Finney et al 2018, thus indicating that there is low confidence for the climate-induced changes of lightning NO x emissions on surface ozone. For wildfires there are a number of uncertainties to be considered as fires modify the structure and distribution of the terrestrial ecosystem with feedbacks in albedo, ecosystem properties, and emissions of BVOCs, GHGs, and other trace species (Voulgarakis and Field 2015, Hantson et al 2016, Li et al 2017, Rabin et al 2017 while there is also feedback of wildfire-induced ozone pollution to global terrestrial gross primary productivity (Yue and Unger 2018). Although regional stagnation and heat waves have been often associated with high ozone concentrations and positive ozone-temperature sensitivity in a number of studies (Fiore et al 2012, Gao et al 2013, Jing et al 2017, Schnell and Prather 2017, Garrido-Perez et al 2018, Zhang et al 2018, weak correlations between observed summertime surface ozone and high-ozone events and the air stagnation index (ASI) or the number of stagnant days over the US have been reported by others (Oswald et al 2015, Sun et al 2017, Kerr and Waugh 2018.
It has been also shown that surface ozone changes, robustly attributable to climate warming alone, may not emerge above the noise driven by internal variability until after 2050 and hence a simulation framework over longer periods is necessary to disentangle climate change from interannual variability (Barnes et al 2016, Lacressonnière et al 2016, Garcia-Menendez et al 2017.
This study provides, for the first time to our knowledge, a global multi-model perspective of climate change penalty and benefit on regional surface ozone based on quantitative estimates from five CMIP6 models of the spatially distributed ∆O 3 /∆T index and the ozone changes for different warming levels. Section 2 presents the data used and the methodology applied in this study. In section 3 the key results of this study are presented and discussed, while, finally, in section 4 the main conclusions are summarized and discussed.

Data and methodology
AerChemMIP (Aerosol Chemistry Model Intercomparison Project), which is endorsed by CMIP6, aims at quantifying the impacts of aerosols and chemically reactive gases on climate and air quality . As part of AerChemMIP, the experiments ssp370SST and ssp370pdSST were carried out The simulations under AerChemMIP ssp370SST and ssp370pdSST experiments were used to quantify the effect of climate change on surface ozone. The ssp370SST experiment follows the ssp370 scenario with time-varying SSTs from coupled model ssp370 simulations. The ssp370pdSST experiment has similar experimental set up as ssp370SST except that sea surface temperatures (SST) and sea ice concentrations are taken from a present-day climatology (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014) from each model's historical simulation and remain constant over the entire period of the experiment (2015-2100). Both ssp370SST and ssp370pdSST experiments consider future anthropogenic emissions following the SSP3_7.0 scenario as described in Turnock et al (2020) and Griffiths et al (2021). Practically, ssp370pdSST simulations depict the effect of emission changes only while ssp370SST simulations depict the combined effect of climate and emission changes. Hence, by subtracting ssp370pdSST from ssp370SST the effect from climate change can be deduced. This is manifested through changes in chemistry, transport, natural emissions, and deposition with the latter two depending on the level of interactiveness at which these processes are represented in the models.
For each model the change in surface ozone concentrations is calculated as the difference between the ssp370SST and ssp370pdSST experiments which can then be averaged to provide the multi-model mean change. The difference between the ssp370SST and ssp370pdSST in global mean annual near-surface temperature was calculated for each model and used as an index of global warming. To assess the changes in surface ozone concentrations for each scenario (ssp370SST and ssp370pdSST), future concentrations from each model have been compared to mean values calculated over the 2005-2014 period from the respective historical AerChemMIP experiment (histSST).
The chemical ozone production and loss terms are provided as model diagnostics (diagnostic variable names: 'o3prod' and 'o3loss') for the calculation of the net chemical ozone production term which can then be related to global warming in the aforementioned model simulations. The ozone production term is based on the sum of all the HO 2 + NO and RO 2 + NO reactions, while the loss term is based on the sum of the following ozone loss reactions; O( 1 D) + H 2 O, O 3 + HO 2 , O 3 + OH and O 3 + alkenes.
Furthermore, the stratospheric ozone tracer diagnostic (diagnostic variable name: 'o3ste') available for two out of the five models (EC-Earth3-AerChem and GFDL-ESM4) was also used in the analysis (see models' description in appendix). In this way, ozone that originates in the stratosphere can be traced through the troposphere. Three out of the five ESMs have interactive BVOC emissions (GISS-E2-1-G, GFDL-ESM4, and UKESM1-0-LL) implying BVOC emission changes under SSP370 global warming. However, BVOC emissions are implemented with a varying level of complexity in the model simulations, as CO 2 inhibition of isoprene emissions is considered only in UKESM1-0-LL. Lightning NO x emissions are calculated interactively in all models. Wildfire emissions are prescribed in the model simulations and hence there is no feedback of wildfire-induced ozone pollution under global warming. Land-use and landcover are also fixed for all models in ssp370SST and ssp370pdSST experiments.
Relevant information for each model and the reference of the simulation experiments are shown in table 1. A short description for each one of the earth system models and chemistry climate models used in this study is provided in appendix.

Results
Based on each model's historical simulation (histSST) the multi-model mean global surface ozone concentration is 32.3 ± 5.8 ppbv over the period 2005-2014. The respective values over the decade 2015-2024 based on ssp370SST and ssp370pdSST experiments are 31.7 ± 4.2 ppbv and 32.0 ± 4.3 ppbv, respectively. The annual cycle of surface ozone (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014) for each model and the multi-model mean is compared with observations from the Tropospheric Ozone Assessment Report across different world regions in figure S1 (available online at stacks.iop.org/ERL/17/024014/mmedia), following that done in Turnock et al (2020). Across most regions, the multi-model mean annual cycle captures the key features of the observed seasonality but with commonly positive biases which become larger in summer than in winter. Especially in oceanic regions, surface ozone is overestimated in the models across all seasons. Additionally, the timing of peak ozone over continental Northern Hemisphere locations occurs earlier in the observations (springtime) than in the models (spring and summer), which is consistent with previous model evaluation studies (Young et al 2018 and references therein). The models' diversity is larger in the continental regions than in oceanic regions with UKESM1-0-LL and EC-Earth3-AerChem diverging the most from the multi-model mean over the Northern Hemisphere. The possible sources of models' structural uncertainty that lead to ozone differences are discussed by Turnock et al (2020) and Griffiths et al (2021). Furthermore, figure S2 shows annual mean surface ozone concentrations and the associated standard deviation for each model based on ssp370pdSST simulations from 2015 to 2100. Despite the diversity among the models on the absolute annual surface ozone concentrations, the models show common patterns in the spatial distribution with higher values simulated in the Northern Hemisphere continental regions and lower over the oceanic regions in the tropics and the Southern Hemisphere. The models also show similar patterns in the spatial distribution of the standard deviation indicating common sources for the simulated ozone variability.
The change in global mean annual near-surface temperature (figure 1) serves as an index of global warming for the SSP3-7.0 scenario with respect to the present day (2005-2014) (ssp370SST-ssp370pdSST) during the time period from 2015 to 2100 for each model and the multi-model mean. Among the models, UKESM1-0-LL presents the highest warming trend over the 21st century, due to its high equilibrium climate sensitivity (5.4 • C based on the method of Gregory et al 2004) as documented by Sellar et al (2019). The smallest warming trend is simulated by the GFDL-ESM4 due to its smaller equilibrium climate sensitivity (2.7 • C based on the method of Gregory et al 2004) as documented by  Dunne et al (2020). The global annual surface ozone change (ssp370SST-ssp370pdSST), shown in figure 1, depicts purely the climate change effect with a prominent decreasing trend over the period from 2015 to 2100 for each model and the multi-model mean.
As it is illustrated in figure 2, there is a prominent anti-correlation between global temperature and global ozone changes for each model and the multimodel mean. The calculated regression coefficients ∆O 3 /∆T for each model are shown in table 2 with a value of −0.96 ± 0.07 ppbv • C −1 for the multi-model mean.
To spatially decompose the regression coefficient ∆O 3 /∆T during the time period from 2015 to 2100, we applied this regression methodology for every model grid point to produce the spatial pattern of the regression coefficient of annual surface ozone change (ssp370SST-ssp370pdSST) over annual surface temperature change (ssp370SST-ssp370pdSST) (figure 3(a)). Similarly, the spatial patterns of the regression coefficient ∆O 3 /∆T were also calculated for the boreal winter/austral summer period including the months December, January, February (DJF) and for the boreal summer/austral winter period including the months June, July, August (JJA), as shown in figures 3(b) and (c), respectively. Figure 3 shows a statistically significant and robust future climate change signal on surface ozone over the oceanic regions ranging between −1 to −2 ppbv • C −1 on an annual basis, as well as for DJF and JJA. Similarly, land-based regions that are located away from large pollution sources also show a comparable (but smaller) reduction in surface ozone due to the climate change signal. However, it should noted the stronger models diversity over the continental areas with non-robust signal. Generally, over regions remote from pollution sources (except the Arctic), there is a consistent decline in mean annual surface ozone concentration due to the global warming associated with ssp370, with the decrease varying spatially from −0.2 to −2 ppbv • C −1 with the strongest decline seen over tropical oceanic regions. This can be attributed to increased net ozone chemical destruction by water vapour under increasing global near-surface temperatures over the 21st century in the SSP3-7.0 scenario. The implication is that over regions remote from pollution sources there is a consistent climate change benefit for baseline ozone due to global warming.
However, over regions close to anthropogenic pollution sources as assumed in the SSP3-7.0 scenario (south-eastern China and India) or close to enhanced natural BVOC emissions sources (e.g. areas in Africa and south America), there are positive values of the regression coefficient ∆O 3 /∆T ranging regionally from 0.2 to 2 ppbv • C −1 implying a regional surface  Table 2. Regression coefficient ∆O3/∆T and explained variance R 2 deduced from linear regression between global annual surface ozone change (ssp370SST-ssp370pdSST) and global annual surface temperature change (ssp370SST-ssp370pdSST) during the time period from 2015 to 2100, for each model and the multi-model mean. ozone penalty due to global warming. The individual models show this robustly for south-eastern China and India (four out five models agree on the sign of change) but there are model differences in regions of Europe and the US as well as in Africa and south America (figure S3 in supplementary material). The magnitude of the above-mentioned regional surface ozone penalty is comparable to the magnitude of the large-scale baseline ozone benefit and so it provides some degree of offsetting on larger continental scale, but regionally the ozone penalty over populated regions would be detrimental. The robust surface ozone penalty due to global warming for south-eastern China and India can be attributed to the much stronger and spatially extensive NO x emissions over these regions in relation to the rest of the world during the 21st century under the SSP370 scenario (figure S4 in supplementary material). The positive effect of climate change on surface ozone over regions close to anthropogenic pollution sources (as assumed in the SSP3-7.0 scenario) is emphasised and explained by the positive local net ozone production rates due to global warming as illustrated in figure S5 (in supplementary material) with the highest robust values over south-eastern China and India, as well as in parts of Africa and south America. This is a robust signal in each one of the models except for EC-Earth3-AerChem in which peak NO x levels in polluted regions are moderated by the relatively coarser resolution (3 • × 2 • ) compared to the other models.

Regression coefficient
Over the US, there is a robust annual mean ozone benefit ranging from −0.2 ppbv • C −1 over the eastern US to −0.8 ppbv • C −1 over the western US which becomes slightly positive but non-robust among the models (an ozone penalty of 0.2 ppbv • C −1 ) over the north-eastern US during JJA (figure 3). Similarly, over Europe, there is a non-robust slight ozone benefit on an annual basis and in DJF which becomes slightly positive but non-robust among the models (an ozone penalty of 0.2 ppbv • C −1 ) during JJA (figure 3). These slightly negative or positive (in JJA) signals over the US and Europe can be attributed: (a) to the generally much lower anticipated NO x emissions (compared to east China and India) during the 21st century under the SSP370 scenario (figure S4 in supplementary material), (b) to contributions from enhanced BVOC emissions for the models that have interactive BVOC emissions (figure S6 in supplementary material), enhanced lightning NO x emissions (figure S7 in supplementary material), and STE (figure S8 in supplementary material), and (c) the competitive role of the large scale background ozone chemical destruction linked to more water vapour under global warming. Hence, these factors (due to compensating effects) lead to a low sensitivity of surface ozone to climate change and a lack of model agreement in the sign of changes under global warming over these regions. Among the models, UKESM1-0-LL reveals the largest ozone penalty over Europe (and North America to a certain extent) in JJA (figure S3 in supplementary material) due to higher NO x emissions over north-western Europe and the north-eastern US as well as due its high equilibrium climate sensitivity. It is interesting to note that the slight ozone penalty values in JJA for Europe and the US are consistent with the positive but non-robust effect of global warming on local net ozone production rate over parts of the US and Europe during JJA (figure S5 in supplementary material).
Over areas of Africa and South America with enhanced natural BVOC emissions sources (figure S6 in supplementary material), there are non-robust signs of positive ozone responses to global warming (figure 3), which are presumably linked to model differences in the parameterization of BVOC emission. As mentioned in section 2, only three out of the five models have interactive BVOC emissions (GISS-E2-1-G, GFDL-ESM4, and UKESM1-0-LL), while among these three models CO 2 inhibition of isoprene emissions is accounted only in UKESM1-0-LL. The models that have interactive BVOC emissions (GISS-E2-1-G, GFDL-ESM4, and UKESM1-0-LL) show increases in BVOC emissions under SSP370 warming over areas in Africa and south America (figure S6 in supplementary material).
Beyond model differences in BVOC emissions, unravelling the effect of global warming on surface ozone over areas of Africa and South America becomes more complicated due to model differences in the contribution of climate change induced lightning NO x emissions (figure S7 in supplementary material). GISS-E2-1-G, MRI-ESM2-0 and UKESM1-0-LL show consistent increases under SSP370 warming especially over areas in Africa and south America while this is not the case for GFDL-ESM4 and EC-Earth3-AerChem. MRI-ESM2-0 (which does not have interactive BVOC emissions) is more sensitive as the positive ozone response is strong and more extensive over areas in Africa and South America (figure S3), which is presumably linked to increases in lightning NO x emissions under SSP3-7.0 warming (figure S7) and more sensitive net chemical ozone production term over land. On the other hand, EC-Earth3-AerChem (which does not have either interactive BVOC emissions) shows weak positive (or even negative) ozone response signal over these regions in Africa and south America (figure S3), while the lightning NO x emission changes under SSP370 warming are also weak.
Aiming to find out the ozone change that would occur due to reaching a warming level such as that targeted for Paris Agreement, we followed an additional approach as illustrated in figure 4. Figure 4 depicts the effect of climate change on surface ozone concentrations for different warming levels using changing global mean surface air temperature in the 21st century. Based on this approach, the surface ozone change for each model corresponds to the same level of warming which is reached at a slightly different time period in each model. This means that there are two levels of inter-model diversityone from the diversity in climate sensitivity causing models to reach warming target at different times and the other is from the sign and magnitude of ozone change simulated for these targets. Similarly, to figure 3(a), this figure illustrates a robust reduction in annual mean surface ozone concentrations in most regions and an ozone penalty over continental regions close to anthropogenic pollution sources (south-eastern China and India) or close to enhanced natural BVOC emissions sources (e.g. areas in Africa and South America, but the changes are not robust) with the signals becoming larger for increasing warming levels. However, the individual CMIP6 models disagree on the sign of the surface ozone response over most of the northern hemisphere area, particularly at the lower surface temperature thresholds, although the agreement between models increases at the 2.5 • C temperature change.
A striking feature in both figures 3 and 4, is the robust sign of increase in surface ozone over the Arctic which is strongest in DJF but is offset by a decreasing signal in JJA to become a weaker positive signal on an annual mean basis. This increase in surface ozone is not accompanied by increase in local net chemical production rate as evident in figure S5 (in supplementary material). However, investigation of the stratospheric ozone tracer diagnostic available for two out of the five models (EC-Earth3-AerChem and GFDL-ESM4) reveals an increase for the stratospheric tracer at near surface over the Arctic for both models in DJF which is also evident with a weaker signal on an annual mean basis (figure S8 in supplementary material). This implies a dynamical cause for the surface ozone increase over the Arctic through stratosphere-to-troposphere transport presumably linked to enhanced STE and stratospheric ozone recovery over the 21st century under global warming.
To compare the climate change effect on surface ozone versus the combined effect of climate and emission changes, the future multi-model change in surface ozone concentrations averaged across the months of boreal summer for ten continental regions and globally in 2050 and 2095 from the five CMIP6 models in the ssp370SST experiments and the difference between the ssp370SST and ssp370pdSST experiments (figure 5). Across most regions, the effect of climate change reduces surface ozone concentrations except for Southern and Eastern Asia, whereas the combined effect of changes in ozone precursor emissions and climate change increases surface ozone concentrations across most world regions. This reveals the dominant role of emission changes on surface ozone changes under future climate change scenarios.

Discussion and conclusions
A multi-model global annual average ozone climate benefit of −0.96 ± 0.07 ppbv • C −1 was calculated solely due to global warming in the ssp370 scenario. Previous modelling studies reported similar global estimates ranging from −0.7 ppbv to −0.88 ppbv for a temperature change of roughly 0.7 • C by 2030s for different scenarios (Dentener et al 2006. In this range, Racherla and Adams (2006) reported a global surface ozone decrease of −1.3 ppbv for a temperature change of 1.7 • C by 2050s. When splitting land and ocean ∆O 3 /∆T, a multi-model mean value of −1.12 ppbv • C −1 for the oceanic regions and −0.44 ppbv • C −1 for the land is estimated. If ∆O 3 /∆T over land is weighted by NO x emissions (shown in figure S4), then a value of +0.06 ppbv • C −1 is calculated, which highlights the importance of NO x emissions over polluted regions in creating an ozone penalty. Nevertheless, the ozone decrease across the oceans dominates the large-scale baseline ozone decrease which through atmospheric circulation provides some degree of offsetting on larger continental scale and may act to moderate any ozone penalty signals over polluted regions.
Spatially, over regions remote from pollution sources, there is a robust decline in mean surface ozone concentrations varying spatially from −0.2 to −2 ppbv • C −1 , with the strongest decline over tropical oceanic regions. A similar conclusion is reached from the inspection of surface ozone concentration changes for different warming levels with the ozone decreases becoming larger with increasing warming levels. The implication is that over regions remote from pollution sources there is a consistent climate benefit for baseline ozone due to global warming. This is mainly linked to the dominating role of enhanced ozone destruction with higher water vapour abundances under a warmer climate, as has been also shown in several previous studies (e.g. In contrast, over regions close to anthropogenic pollution sources (as assumed in the SSP3-7.0 scenario) or close to enhanced natural BVOC emission sources, there are ozone increases with a rate ranging regionally from 0.2 to 2 ppbv • C −1 implying a regional surface ozone penalty due to global warming. The role of local photochemical ozone production for the ozone increase over these regions is justified by positive local net ozone production rates due to global warming. The individual CMIP6 models show this chemically driven ozone penalty robustly for southeastern China and India due to the stronger and spatially extensive NO x emissions over these regions than over other continental regions during the 21st century under the SSP3-7.0 scenario. Ozone penalties at polluted regions (with large NO x emissions, including large biomass burning emissions) have been reported in various previous studies (Hauglustaine et  Over Europe and the US, there are slightly negative or positive ozone changes, which are related to the much lower NO x emissions than in eastern China and India, as well as to contributions from enhanced BVOC and lightning NO x emissions, STE, and the competitive role of the large-scale baseline ozone chemical destruction. Hence, compensating effects from these processes and model differences in parameterizations of BVOC and lightning NO x emissions, and STE lead to low sensitivity of surface ozone to climate change and lack of model agreement for the sign of changes under global warming. The above mentioned heterogeneity in model results is also found in previous regional studies focusing on North America (Gonzalez-Abraham et al 2015, Val Martin et al 2015, Schnell et al 2016, He et al 2018, Nolte et al 2018, Rieder et al 2018 or Europe (Katragkou et al 2011, Colette et al 2015, Lacressonnière et al 2016, Schnell et al 2016, Fortems-Cheiney et al 2017. For European land areas, previous studies have estimated a summertime average ozone climate penalty as a function of the average European surface temperature anomaly of +0.17 ppbv • C −1 (Colette et al 2015), which is consistent with the estimated range of 0 to +0.2 ppbv • C −1 inferred by the multi-model mean of CMIP6 models over parts of Europe. Furthermore, positive relationships were also reported for three CCMs, with the climate penalty ranging between 0.34 and 1.20 ppbv • C −1 for three midlatitude regions in North America, Europe and East Asia from May to September (Doherty et al 2013).
The low sensitivity of surface ozone to climate change over Europe and the US versus the high sensitivity over south-eastern China and India points to the fact that climate change enhances the efficiency of precursor emissions to generate surface ozone in polluted regions, but the magnitude of this effect depends on the regional emissions considered versus the competitive role of the large-scale baseline ozone chemical destruction. Most of previous studies quantified ozone climate penalty for different global warming scenarios while holding the ozone precursor emissions fixed at present day (Colette et al 2015, Schnell et al 2016. Schnell et al (2016) showed climate-driven increases in summertime surface ozone over the north-eastern US, parts of Europe, and north-eastern China under RCP8.5 warming scenario with anthropogenic precursor emissions fixed to present day. In our study, anthropogenic precursor emissions follow the business-as-usual scenario SSP3-7.0. However, the use of either present-day or future business-as-usual scenarios of anthropogenic precursor emissions could bias the climate ozone penalty results. This has been pointed out by Colette et al (2013) reporting that the climate ozone penalty will decrease in magnitude and even become a net benefit under mitigation scenarios for anthropogenic precursor emissions. Hence it is important that future studies quantifying climate ozone penalty should consider both business-as-usual and mitigation scenarios for anthropogenic precursor emissions. Furthermore, the amplitude of the climate change penalty on ozone over polluted regions may be different in high-resolution (regional and urban scale) models in comparison to coarse resolution global models as several controlling processes are resolution dependent, including e.g. meteorological conditions, local emissions, and sensitivity to the chemical regime (VOC limited versus NO x limited) (Lauwaet et al 2014, Markakis et al 2014, 2016.
Ozone increases due to global warming are also revealed in regions close to enhanced natural BVOC emission sources such as in Africa and south America, but the signal is not robust in the multi-model mean as only three out of the five models have interactive BVOC emissions while there also model differences in the parameterization of BVOC emission. Nevertheless, these three models (with interactive BVOC emissions) show increases due to global warming over these regions. Previous studies have also shown increase in BVOC emissions because of climate warming on local to global scales (Weaver et al 2009, Jiang et al 2018. However, a growing number of recent studies indicates that the CO 2 inhibition of biogenic isoprene emissions reduces future increases in BVOC emissions and might even result in a decrease (Tai et al 2013, Fu and Liao 2016, Hollaway et al 2017, Lin et al 2020. In our study CO 2 inhibition of isoprene emissions is accounted only in UKESM1-0-LL simulations. Furthermore, there are feedback processes related to biosphereatmosphere interactions which are not considered in the presented simulations. For example, ozone itself interacts with vegetation either via stomatal uptake with negative feedback in BVOC emissions and ozone or through a reduced deposition velocity of ozone, resulting in positive feedback on ozone (Franz et al 2017, Sadiq et al 2017, Lombardozzi et al 2018, Zhou et al 2018. Other relevant studies indicated that the dry deposition of ozone, driven by non-stomatal mechanisms (rather than by stomatal uptake) may alleviate the positive feedback on ozone over many vegetated surfaces (Clifton et al 2014, Horváth et al 2018. Even if the emissions of biogenic VOCs do increase with climate warming, the associated response of surface ozone over BVOC-rich areas is still quantitatively uncertain due to the different model assumptions on the yields of isoprene nitrates and their subsequent NO x -recycling ratios (Squire et al 2015). Furthermore, there is lack of knowledge on quantifying the effect of natural methane emissions on surface ozone under future climate change. The effect of climate change on ozone via changes in methane concentrations and lifetime are not considered as the models use the same prescribed methane concentrations in all simulations. Overall, the integrated effect of biosphere interactions on ozone remains poorly understood and there are model discrepancies for their impact on surface ozone under global warming with even opposite signs of ozone changes (Squire et al 2015, Val Martin et al 2015, Schnell et al 2016, Pommier et al 2018. There is a robust sign of increase in surface ozone over the Arctic in DJF by roughly 1 to 2 ppbv • C −1 which seems to be associated with dynamical causes through STE processes under global warming as the stratospheric ozone tracer (in the two available models) at the surface reveals an average increase of about 1 to 2 ppbv over the whole period 2015-2100. A significant increase of surface ozone due to climate change over the Arctic has been also reported by Zeng et al (2008) in February/March. The changes of the net stratospheric influx in STE are linked to changes of the stratospheric Brewer-Dobson Circulation and the amount of ozone in the lowermost stratosphere, which are strongly influenced in a changing climate by the emissions of ODSs and GHGs (Butchart 2014, Morgenstern et al 2018. Total column ozone (reflecting mostly stratospheric ozone) is projected to return to 1960s values by the middle of the 21st century under the SSP2-4.5, SSP3-7.0, SSP4-3.4, SSP4-6.0 and SSP5-8.5 scenarios (Keeble et al 2021). Furthermore, analysis of five CMIP6 models (UKESM1-0-LL, CESM2-WACCM, GFDL-ESM4, MRI-ESM2-0, and GISS-E2-1-G) showed pronounced stratospheric ozone increases, which impact the tropospheric-ozone abundance through STE in the SSP3-7.0 scenario across all models (Griffiths et al 2021). Our findings are in line with a number of previous model studies showing that both enhanced stratospheric ozone influx into the troposphere and stratospheric ozone recovery will tend to increase the future tropospheric ozone levels, affecting surface ozone by rougly 1 and 2 ppbv in the extratropics (Sekiya and Sudo 2014, Hess et al 2015, Banerjee et al 2016, Meul et al 2018, Akritidis et al 2019. However there are still uncertainties and a large spread among the different models for the effect of STE changes on tropospheric and surface ozone related to stratosphere-troposphere coupling and tropospheric mixing processes among the models (Morgenstern et al 2018).
Apart from the limitations resulting from the uncertainties in processes such as STE and atmosphere-biosphere interactions discussed above, we should also consider the limitations in quantifying the effects on surface ozone under a changing climate from lightning NO x emissions as well as from wildfires. Lightning NO x emissions are calculated interactively in the model simulations of this study, but there are model differences in climate change induced lightning NO x emissions (ssp370SST-ssp370pdSST differences). Furthermore, the direction of change of lightning activity in the future remains also highly uncertain (Clark et al 2017, Finney et al 2018. In our study, wildfire emissions are not interactive in the model simulations and hence there is no feedback of wildfire-induced ozone pollution under global warming. Another source of uncertainty of our study is driven by the diversity in models' equilibrium climate sensitivity which leads to diversity in reaching different warming levels and the associated surface ozone changes. The comparison of the climate change impact effect on surface ozone versus the combined effect of climate and emission changes reveals the dominant role of emission changes controlling surface ozone changes under future climate change. This is in line with the assessments of IPCC AR5 and AR6, which concluded with high confidence that the response of surface ozone to climate-driven changes is more uncertain than the response to emission-driven changes (Kirtman et al 2013, Szopa et al 2021. In AR6, the evolution of surface ozone is shown to be very variable within each region depending on the choices assumed in the various SSP scenarios regarding air pollution control (Szopa et al 2021).

Key-note messages
Overall, a warmer climate under SSP3-7.0 scenario is expected to reduce surface ozone (ozone benefit) in unpolluted regions (from −0.2 to −2 ppbv • C −1 ) and globally (approximately −1 ppbv • C −1 ) because of greater water vapour abundance accelerating ozone chemical loss. It is also expected to increase surface ozone (ozone penalty) by a few ppbv (up to 2 ppbv • C −1 ) over polluted regions depending on the regional emissions considered in this study within the SSP3-7.0 scenario and the competitive role of the large-scale baseline ozone chemical destruction. However, there are uncertainties in several processes affected in a warmer climate including biosphereatmosphere interactions, STE, wildfires and lighting NO x emissions, which can affect and modify future baseline and regional/local surface ozone levels. The overall response of surface ozone to climate-driven changes is more uncertain than the response to emission-driven changes, which play the dominant role in controlling surface ozone changes under future climate change. It is important for air pollution policies that upcoming studies quantifying climate ozone penalty should consider both business-asusual and different mitigation scenarios for anthropogenic precursor emissions. This work underlines an additional need to reduce ozone pollution to avoid or limit the climate penalty at polluted regions due to a warmer climate.
Our multi-model analysis focuses at large-scale responses of annual and seasonal surface ozone to climate change on decadal climatic timescales. Thus, due to the non-linear nature of ozone chemistry the investigation of climate-driven changes in pollution episodes and ozone extremes at hourly and daily temporal resolution is masked. Further work on ozone extremes in selected polluted regions based on ESMs would be an asset for air pollution policies as heat waves and air pollution episodes pose a serious threat to human health and may worsen under future climate change depending on the scenario.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: https://esgf-node.llnl.gov/search/cmip6/.

Author contributions
P Z is the main author drafting the paper and contributed to the interpretation of data analysis. D A and S T contributed to data analysis, figure preparation and interpretation of the results. S S and V N motivated the analysis design of the paper and contributed to the interpretation of data analysis. A K G contributed to data analysis and figure preparation. S E B and K T contributed in the GISS-E2-1-G simulations; M D and N O contributed in the MRI-ESM2-0 simulations; T v N and P L S contributed in the EC-Earth3-AerChem simulations; J K and F M O contributed in the UKESM1-0-LL simulations; L W H contributed in the GFDL-ESM4 simulations. All authors commented on the paper. We would also to thank the reviewers for their constructive comments.

A.1 GFDL-ESM4
The atmospheric component of GFDL-ESM4 (Dunne et al 2020), called AM4.1, includes an interactive tropospheric and stratospheric gas-phase and aerosol chemistry scheme . The model includes 76 prognostic (transported) chemical tracers and 40 diagnostic (non-transported) chemical species, with 43 photolysis reactions, 190 gasphase kinetic reactions, and 15 heterogeneous reactions. The tropospheric chemistry includes reactions for the NO x -HO x -O x -CO-CH 4 system and oxidation schemes for other NMVOCs. The stratospheric chemistry accounts for the major ozone loss cycles (O x , HO x , NO x , ClO x , and BrO x ) and heterogeneous reactions on liquid and solid stratospheric aerosols as in Austin et al (2012). The chemical system is solved using an implicit Euler backward method with Newton-Raphson iteration. Photolysis rates are calculated interactively using the FAST-JX version 7.1 code, accounting for the radiative effects of simulated aerosols and clouds. Emissions of BVOCs, including isoprene and monoterpenes, are calculated online in AM4.1 using the Model of Emissions of Gases and Aerosols from Nature (Guenther et al 2006), as a function of simulated air temperature and shortwave radiative fluxes. Details on the chemical mechanism are included in Horowitz et al (2020).
Lightning NO x emissions are calculated interactively as a function of subgrid convection diagnosed by the double-plume convection scheme described by Zhao et al (2018). The lightning NO x source is calculated as a function of convective cloud-top height, following the parameterization of Price et al (1997) and is injected with the vertical distribution of Pickering et al (1998). Anthropogenic and biomass burning emissions are prescribed from the dataset of Hoesly et al (2018) and van Marle et al (2017) developed in support of CMIP6. Natural emissions of ozone precursors are not calculated interactively but they are prescribed as in Naik et al (2013). The stratospheric ozone tracer diagnostic (o3ste) is defined using an e90 tropopause (Prather et al 2011). Above the local e90 tropopause, o3ste is set to be equal to ozone for each time step. Below the e90 tropopause, o3ste is lost with a loss frequency corresponding to odd-oxygen loss via reactions O 1 D + H 2 O, HO 2 + O 3 , OH + O 3 , O 3 + alkenes, NO 2 + OH, heterogeneous loss of N 2 O 5 , NO 3 , and NO 2 , and dry deposition.
A.2 UKESM1-0-LL UKESM1-0-LL is the UK's Earth System Model (Sellar et al 2019). It is based on the Global Coupled 3.1 (GC3.1) configuration of HadGEM3 (Williams et al 2018), to which various Earth system components have been added e.g. ocean biogeochemistry, terrestrial carbon/nitrogen cycle, and atmospheric chemistry. The atmospheric and land components are described in Walters et al (2019). The chemistry scheme included is a combined stratospheretroposphere chemistry scheme (Archibald et al 2020) from the UK Chemistry and Aerosol (UKCA) model, combining the stratospheric chemistry scheme of Morgenstern et al (2009) with the tropospheric chemistry scheme of O'Connor et al (2014). The horizontal resolution is 1.25 × 1.875 (∼140 km at mid-lats) with 85 levels in the vertical (up to 85 km) and 84 chemical tracers used to simulate chemical cycles of O x , HO x , and NO x , as well as oxidation reactions of CO, CH 4 , and NMVOCs. CH 4 is prescribed by surface concentrations from each CMIP6 scenario. Anthropogenic and biomass burning emissions are prescribed (van Marle et al 2017, Hoesly et al 2018), but emissions of isoprene and monoterpenes are interactive based on the interactive biogenic VOC emission model (Pacifico et al 2011). Lightning emissions of NO x (LNO x ) are also interactive using the cloud top height parameterization of Price and Rind (1992). Other natural emissions are prescribed as climatologies as also discussed fully in Archibald et al (2020). For volcanic eruptions, internally consistent stratospheric AODs and SADs are prescribed for both the volcanic forcing and for the UKCA stratospheric heterogeneous chemistry. The UKESM1-0-LL set up for the ssp370pdSST was to fix sea surface temperatures, sea ice fields, CO 2 concentrations, and surface water dimethyl sulfide and chlorophyll concentrations at present day values. Fixing CO 2 in ssp370pdSST would not allow the act of CO 2 inhibition in the absence of a temperature change so that BVOC emissions would stay close to present-day levels. MRI-ESM2-0 includes interactive chemistry and aerosols in the atmosphere. The atmospheric chemistry model, MRI-CCM2.1, calculates evolution and distribution of the ozone and other trace gases in the troposphere and middle atmosphere Shibata 2011, Yukimoto et al 2019a). The model includes 64 prognostic chemical species and 24 diagnostic chemical species, with 184 gas-phase reactions, 59 photolysis reactions, and 16 heterogeneous reactions. It considers O x -HO x -NO x -CH 4 -CO chemical system and NMVOC oxidation reactions, as well as the major stratospheric chemical system. Anthropogenic and biomass burning emissions are prescribed (van Marle et al 2017, Hoesly et al 2018. Lightning emissions of NO x are diagnosed at 6 h intervals following the parameterization of Price and Rind (1992). Other natural emissions such as biogenic, soil, and ocean emissions are prescribed as climatologies (Deushi and Shibata 2011).

A.4 EC-Earth3-AerChem
EC-Earth3-AerChem is essentially EC-Earth3 extended with an additional component to simulate aerosols and atmospheric chemistry (Döscher et al 2021, van Noije et al 2021. The atmospheric GCM is based on IFS cycle 36r4, which includes the land surface model HTESSEL. Its horizontal resolution is TL255 (triangular truncation at wavenumber 255 in spectral space with a linear N128 reduced Gaussian grid, corresponding to a spacing of about 80 km). Its atmospheric grid consists of 91 layers in the vertical direction and has a model top at 0.01 hPa. The time step applied in IFS is 45 min. For EC-Earth3-AerChem, the CMIP6 forcings prescribed in IFS are the solar forcing (Matthes et al 2017), well-mixed greenhouse gas concentrations (CO 2 , N 2 O, CFC-12, and CFC-11 equivalents; Meinshausen et al 2017Meinshausen et al , 2020, and stratospheric aerosol radiative properties. The evolution of methane is constrained by the surface mixing ratio data provided by CMIP6 (Meinshausen et al 2017(Meinshausen et al , 2020. Aerosols and atmospheric chemistry are simulated with the Tracer Model version 5 (TM5), specifically release 3.0 of the massively parallel version of TM5 (TM5-mp 3.0). It runs at a horizontal resolution of 3 • × 2 • (longitude × latitude) with 34 layers in the vertical direction. TM5 was integrated as a module coupled to IFS within EC-Earth (van Noije et al , 2021. The chemistry scheme of TM5 accounts for gas-phase, aqueous-phase and heterogeneous chemistry (van Noije et al 2021) with the gas-phase reaction scheme being a modified version of the CB05 carbon bond mechanism. The mixing ratios of ozone in the stratosphere are nudged towards zonal mean fields calculated from the three-dimensional input data sets provided by CMIP6 ( (Gidden et al 2019). The amount of NO x produced in lightning discharges is calculated interactively whereas other natural emissions of ozone precursors are prescribed, as documented in van Noije et al (2021). The stratospheric ozone tracer diagnostic (o3ste) above the model level close to 140 hPa is set to be equal to ozone for each time step. In the troposphere, o3ste has no tropospheric chemical production but its loss through O( 1 D) + H 2 O, HO 2 + O 3 , OH + O 3 and dry deposition is considered.
A.5 GISS-E2-1-G GISS-E2-1-G is the NASA Goddard Institute for Space Studies (GISS) chemistry-climate model version E2.1 using the GISS210 Ocean v1 (G01) model. The model configurations submitted for CMIP6 are described in detail by Kelley et al (2020) and Miller et al (2021). The atmospheric component was run with horizontal resolution of 2.5 • × 2 • (longitude × latitude) with 40 hybrid sigma-pressure vertical layers extended from the surface to 0.1 hPa. Online interactive chemistry follows the GISS Physical Understanding of Composition-Climate INteractions and Impacts (G-PUCCINI) mechanism for gas-phase chemistry (Shindell et al 2013 and either the One-Moment Aerosol (OMA) or the Multiconfiguration Aerosol TRacker of mIXing state (MATRIX) model for the condensed phase (Bauer et al 2020). The gas-phase mechanism includes 146 reactions (including 28 photodissociation reactions) acting on 47 species throughout the troposphere and stratosphere including five heterogeneous reactions. The model advects 26 (OMA) or 51 (MATRIX) aerosol particle tracers and 34 gasphase tracers. Anthropogenic and biomass burning emissions are prescribed following the CMIP guidelines. NO x emissions are calculated online in deep convection as described by Kelley et al (2020). Soil microbial NO x emissions are prescribed from climatology. Biogenic emissions of isoprene are calculated online and respond to temperature , but are prescribed for alkenes, paraffins and terpenes. Methane is prescribed as a surface boundary condition but allowed to advect and react with the chemistry in the historical runs and a subset of the SSP simulations; some future simulations used interactive online methane emissions following Shindell et al (2004).