Decadal O 3 variability at the Mt. Cimone WMO/GAW global station (2,165 m a.s.l., Italy) and comparison with two high-mountain "reference" sites in Europe

Tropospheric ozone (O 3 ) is a greenhouse gas as well as a harmful air pollutant with adverse effects on human health and vegetation: The observation and attribution of its long-term variability are key activities to monitor the effectiveness of pollution reduction protocols. In this work, we present the analysis of multi-annual near-surface O 3 (1996–2 0 16) at the Mt. Cimone (CMN, Italian northern Apennines) WMO/GAW global station and the comparison with two “reference” high-mountain sites in Europe: Jungfraujoch (JFJ, Swiss Alps) and Mt. Zugspitze (ZUG/ZSF, German Alps). Negative O 3 trends were observed at CMN over the period 1996–2 0 16 (from – 0 .19 to – 0 .22 ppb yr –1 ), with the strongest tendencies as being observed for the warm months (May–September: – 0 .32 ppb yr –1 during daytime). The magnitude of the calculated O 3 trends at CMN are 2 times higher than those calculated for ZUG/ZSF and 3–4 times higher than for JFJ. With respect to JFJ and ZUG/ZSF, higher O 3 values were observed at CMN during 2 00 4–2 00 8, while good agreement is found for the remaining periods. We used Lagrangian simulations by the FLEXPART particle dispersion model and near-surface O 3 data over different European regions, for investigating the possibility that the appearance of the O 3 anomalies at CMN could be related to variability in the atmospheric transport or in near-surface O 3 over specific source regions. Even if it was not possible to achieve a general robust explanation for the occurrence of the high O 3 values at CMN during 2 00 4–2 00 8, the variability of (1) regional and long-range atmospheric transport at CMN and (2) European near-surface O 3 could motivate the observed anomalies in specific seasons and years. Interestingly, we found a long-term variability in air mass transport at JFJ with enhanced (decreased) contributions from Western European (intercontinental regions).


Introduction
The Mediterranean basin is recognized as a global hot spot region for climate change (Giorgi and Lionello, 2008) and air quality (Monks et al., 2009;Stocker et al., 2013). During summer, this region is exposed to high levels of tropospheric ozone (O 3 ), which can have adverse effects on population health and ecosystems also playing an impact on regional climate as short-lived climate forcer (United Nations Environment Programme and World Meteorological Organization, 2011). Due to its absorbing properties in the long-wave radiation, tropospheric O 3 is a fundamental component in the Earth's climate. Moreover, it acts as a powerful oxidizing agent able to determine the overall oxidative capacity of the troposphere and determining the fate of other atmospheric compounds. Due to its high oxidizing capacity, it can cause serious health problems, especially respiratory illnesses and cardiovascular diseases, leading to premature death in some cases (Fleming et al., 2018, and reference therein). It also damages vegetation such as forests and agricultural crops and the services they provide, in particular production and carbon sequestration (Stocker et al., 2013).
Quantitative understanding of long-term baseline O 3 variability can provide valuable information for a number of tasks, for example, in support of regional strategies for the control of tropospheric O 3 as well as to contribute to the validation of global and regional air quality or climate models. Indeed, process-oriented assessments of atmospheric models are needed to build confidence in their utility for assessing the effectiveness of pollution control strategies or projecting pollution extremes under future climate scenarios (e.g., Monks et al., 2015).
The central part of the Mediterranean basin is characterized by the presence of large urban areas and distributed source regions of O 3 precursors (e.g., the Po basin). As shown by satellite measurements, the central Mediterranean basin (i.e., from 5 E to 20 E) is the region where O 3 is maximized in summer at the surface, with the influence of the "eastern" basin free tropospheric O 3 pool still detectable at 1 and 2 km a.s.l. (Safieddine et al., 2014;Zanis et al., 2014).
Since the lifetime of tropospheric O 3 spans from days to weeks and the lifetime of its precursors have an even larger range, air mass transport can significantly affect O 3 variability over specific regions or locations. For example, Cui et al. (2011), by analyzing 19 years of O 3 variability at the Jungfraujoch observatory (JFJ, Switzerland), found a possible impact of reduced NO emissions within the European planetary boundary layer (PBL) on winter and summer O 3 levels.
Climatic variability of atmospheric circulation from the interannual to the decadal scale affects tropospheric O 3 (Lin et al., 2014). For instance, the North Atlantic Oscillation (NAO) has been shown to influence winter and summer climate over the Mediterranean basin (Wang et al., 2011), and due to the effect on circulation patterns and meteorological conditions in the Mediterranean basin, a potential influence of NAO variability on tropospheric O 3 has been suggested (see, e.g., Pausata et al., 2012;Cuevas et al., 2013). Doche et al. (2014) reported that the relative position and the strength of the meteorological systems (Azores Anticyclone and Middle Eastern Depression) over the Mediterranean are key factors in explaining both the variability and the anomalies of O 3 in the lower troposphere of this region. Kalabokas et al. (2017) pointed out that spring high O 3 episodes over the western Mediterranean and central Europe can be linked to specific synoptic meteorological conditions, with the horizontal advection and the subsidence of tropospheric air masses influencing O 3 across the PBL.
Recently, Cooper et al. (2020) analyzed surface O 3 trends at 27 globally distributed remote locations. Since 1995, the Northern Hemisphere sites are nearly evenly split between positive and negative O 3 trends. In particular, negative O 3 trends (-1.8 ppb/decade and -2.1 ppb/decade) were reported for the two high elevation Alpine sites located at Sonnblick (Austria) and Zugspitze (Germany), while a weak positive trend (0.5 ppb/ decade) was reported for JFJ. Parrish et al. (2020) combined the time series from these Alpine sites to one single data set: The analysis of the long-term O 3 changes for this combined data set suggested that the decades-long baseline increase has ended around 2005 over the Alpine region. It should be noted that Cooper et al. (2020) demonstrated that, especially in summer, these Alpine sites frequently sample polluted air from the European boundary layer air. Thus, they can only reflect the lower free tropospheric O 3 , and their longterm trends can be influenced by the signal of the European PBL air masses. Indeed, despite what was observed at the European mountain sites, by exploiting the IAGOS (In-Service Aircraft for a Global Observing System) database (Petzold et al., 2015), Gaudel et al. (2020) reported an increase in median O 3 over Europe from 1 to 3 ppb/decade from the lower to the free troposphere between 1994 and 2016.
Despite the good availability of surface O 3 observations in Europe compared to other regions in the world (Schultz et al., 2017), only very few atmospheric research observatories in the central Mediterranean basin provide continuous high-quality information on baseline "near-surface" O 3 . The longest available time series in the Mediterranean is available at the Global Station "Mt. Cimone," located at 2,165 m a.s.l. in the Italian northern Apennines. In the framework of the Global Atmosphere Watch (GAW) program by the World Meteorological Organization (WMO), near-surface O 3 measurements were continuously carried out since 1996. Previous investigations by Cristofanelli et al. (2015) found no significant O 3 trends for the period 1996-2011, but a slowing down of the positive longterm trend was observed for the summer season. However, larger positive O 3 anomalies were observed for the years 2004-2008 with respect to other baseline measurement sites in the Alps and in the Mediterranean basin. These high O 3 values were, at least tentatively, attributed to the synergic occurrence of the atmospheric processes (i.e., occurrence of heat-waves, deep stratospheric intrusions, and NAO variability).
The goal of this article is to compare the long-term O 3 variability and trends/tendencies at CMN with those observed at two other mountain "reference" stations in Europe as well as to analyze the differences observed during the period 2004-2008 as a function of (1) the O 3 observations carried out at Febbio (FEB), a mountain site located midway between CMN and the Po basin; (2) the changes in the atmospheric transport regimes as deduced by the application of a Lagrangian particle dispersion model; and (3) the interannual variability of near-surface O 3 over Europe as depicted by the Air-Base database. In particular, CMN data series are compared with those observed at the WMO/GAW global stations JFJ (Swiss Alps) and Zugspitze (merged data set from "Gipfel" and "Schneefernerhaus" stations; ZUG/ ZSF, German Alps). In the following, we will use the code "ZUG/ZSF" when referring to the merged data set or to information/studies representative for both the stations, while the specific codes "ZUG" and "ZSF" will be used for referring to the station "Gipfel" (summit) and "Schneefernerhaus," respectively. We decided to use JFJ and ZUG/ZSF as "reference" stations for the O 3 variability on the basis of their spatial proximity to CMN, the availability of multi-decadal O 3 data record and the high level of maturity characterizing their O 3 data sets (see Schultz et al., 2017;Cooper et al., 2020;Parrish et al., 2020, and references therein). A comparison with the IAGOS data set, albeit of interest to evaluate the representativeness of CMN observations with respect to the free tropospheric O 3 levels, is outside the scope of this article and thus is not provided here.

Observation sites and near-surface O 3 measurements
The geographical locations of the considered measurement sites are reported in Figure SM1, while specific details about CMN, JFJ, and ZUG/ZSF O 3 observations and related quality control methodologies can be found elsewhere (e.g., Gilge et al., 2010;Cristofanelli et al., 2015). All three stations have been audited by the WMO/GAW World Calibration Centre for Surface Ozone (WCC-Empa) for complete assessments of QA/QC procedures, measurement systems, and related protocols (Herzog et al., 1996(Herzog et al., , 1997(Herzog et al., , 1999Zellweger et al., 2001Zellweger et al., , 2006aZellweger et al., , 2006bZellweger et al., , 2011Zellweger et al., , 2012Zellweger et al., , 2015Zellweger et al., , 2018 Mt. Cimone (44, 17 N, 10, 68 E, 2,165 m a.s.l.) is the highest peak of the Italian northern Apennines. The nearsurface O 3 measurements presented in this work have been carried out at the "O. Vittori" atmospheric observatory, which is operated by the Institute of Atmospheric Sciences and Climate of the National Research Council of Italy, and it is part of the WMO/GAW Global Station (GAW ID: CMN). As reported by previous investigations, the atmospheric observations carried out at CMN can be considered representative of the free tropospheric conditions of the Mediterranean basin/South Europe during the cold months (see, e.g., Bonasoni et al., 2000), while during the warm season, vertical transport of air masses from the regional PBL is detectable due to the higher vertical extent of the PBL mixing and the activation of thermal wind circulation along the mountain slopes and the valleys (e.g., Cristofanelli et al., 2007). At CMN, O 3 measurements have been carried out by means of UVabsorption analyzers: Dasibi 1108 W/ GEN (1996GEN ( -2012 and Thermo Scientific Tei49i (2012. Starting from 2006, the traceability of O 3 measurements to primary standards has been assured by comparing the laboratory calibrators (Dasibi 1008 PC and, since 2012, Thermo Scientific Tei49i-PS) to reference photometers hosted at the National Institute of Metrological Research (IMGC-O3SRP primary standard) and at the WMO/GAW "World Calibration Centre for surface O 3 , CO, CH 4 , and CO" (WCC-Empa Standard Reference Photometer SRP#15). Previously (from 1996 to 2006), calibrations were directly executed by the instrument manufacturer.

FEB
For interpreting the O 3 variability observed at CMN, observations from an air quality station at FEB (44, 30 N, 10, 43 E; 1,121 m a.s.l.), located 60 km to CMN, are also considered. FEB is a small mountain village (172 inhabitants) located in a rural area. The air quality station managed by ARPAE Emilia-Romagna is located at the border of a forestry area outside the village on the northeastern ridge of Mt. Cusna (2,150 m a.s.l.). No major traffic routes or industries are present near the station, while the next major industry cluster is located 35 km northeast from FEB in the direction of the Po basin: Thus, the station is not influenced by direct emission of specific sources like traffic, industries, or domestic heating. Being located in the Apennines at a lower altitude than CMN, this station is used to track the role of diurnal thermal transport of polluted air masses from the Po basin to the Apennines foothills. Here, we adopt an approach similar to Gheusi et al. (2011) to investigate the diurnal transport of air masses from the boundary layer across a mountain region. Even if not located along the same valley (besides FEB, no other monitoring sites are present at low altitude in the foothills of Apennines), we expect that the O 3 observations carried out at this station can provide useful hints to better understand the diurnal variability of O 3 at CMN with a special emphasis on the role of air mass transport from the Po basin.

JFJ
The high alpine station JFJ (7.98 N; 46.55 E; 3,580 m a.s.l.) is also a global station within WMO/GAW (GAW ID: JFJ). It is located in the northern part of the Swiss Alps and belongs to the first topographical barrier for the frequent westerly winds in central Europe. Its location is relatively remote, with the nearest villages more than 8 km in horizontal and 2.5 km in vertical distance and is only weakly influenced by local anthropogenic sources. A previous assessment of the spatial representativeness of in situ measurement sites in Europe (Henne et al., 2010) revealed that only Lampedusa (central Mediterranean Sea) and Mace Head (Western Ireland) can be considered less influenced by European boundary layer emission than JFJ. Cui et al. (2011) concluded that long-term O 3 changes at JFJ are most likely caused by processes not covered by the 20-day backward trajectories applied, while Logan et al. (2012) reported a good spatial representativeness of the JFJ O 3 data: They showed a good agreement (in particular since the middle of 1990s) with other European mountain sites and data from regular ozone soundings. Due to its height and, thus, its large footprint, JFJ can occasionally be influenced by emissions from a wide area surrounding the Alps. The emission source regions that have a detectable impact on the measurements at JFJ have been investigated, for example, by Reimann et al. (2004) for longlived species and Cui et al. (2011) for O 3 . Reimann et al. (2004) concluded that emissions from an integrated area in central Europe including Switzerland, northern Italy, France, southern and western Germany, the Benelux countries, and to limited extent northeastern parts of Spain can be observed at JFJ. A recent analysis by Cooper et al. (2020) showed that approximately 18% of the air at JFJ is from PBL during winter (December-February), while in summer, this influence increases to 33% during night and 45% during day. Thus, an influence from the European atmospheric boundary layer cannot be completely neglected, and an influence from surface emissions or deposition is expected at JFJ (see also Griffiths et al., 2014).
At JFJ, measurements were carried out by an Environics S300 in 1996, Tei49C (Thermo Environmental Instruments Inc.) instruments from January 1997 to June 2012, and Thermo Scientific Tei49i instruments since June 2012. For most of the time, two instruments were operated in parallel. Instruments were usually replaced in 2 years and were sent to manufacturer representatives in Switzerland for service such as cleaning and replacement of wear parts. The instruments undergo comprehensive tests and quality checks after first receipt and services. At JFJ, instruments were calibrated every 3 months with Tei49C-PS instruments that were traceable by the WCC-Empa SRP#15.
2.1.4. Zugspitze (ZUG/ZSF) Zugspitze (47.4 N,11.0 E) is a WMO/GAW station located in northern flank of the German Alps. It is located approximately 90 km southwest of Munich, while the nearest (10 km far) major town is Garmisch-Partenkirchen (about 27,000 inhabitants, 708 m a.s.l.). Together with CMN, ZUG was identified by Henne et al. (2010) as weakly influenced by surface fluxes from the European PBL, but the works by Yuan et al. (2019) and Carnuth et al. (2002) pointed out a not negligible impact of PBL transport to the atmospheric composition observations. O 3 measurements were carried out by the Fraunhofer Institute for Atmospheric Environment Research (IMK-IFU) at the summit of Mt. Zugspitze (2,962 m a.s.l., GAW ID: ZUG) from 1978 to 2010 (see, e.g., Scheel et al., 1997), while the German Federal Environmental Agency (UBA) established O 3 measurements since 2001 at the lower elevation Schneefernerhaus station (2,671 m a.s.l., GAW ID: ZSF), right between the summit and the skiing area (see, e.g., Gilge et al., 2010). At ZUG, O 3 measurements were carried out by a UV-absorption analyzer (Thermo Environmental Tei49C) since March 1996: For 3 years, parallel measurements were carried out also by the previous chemiluminescence instrument (Bendix 8002), reporting similar performance in terms of accuracy . Also, UBA adopted a UV technique for measurements at ZSF. At the station, different instruments were used during the period 2000-2016 (Thermo Environmental Tei49C, Thermo Scientific Tei49i, and Horiba APOA-370) as well as different O 3 standard (Thermo Environmental Tei49C-PS and Thermo Scientific Tei49i-PS), which used to be calibrated against the UBA standard (SRP#29) on a yearly basis. To compare the Zugspitze O 3 record with CMN and JFJ, we merged the time series at ZUG (1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010) and ZSF (2011ZSF ( -2016 available at the GAW World Data Center for Reactive Gases (WDCRG) hosted at http:// ebas.nilu.no (please note that 2011 data, still not available from WDCRG, have been directly provided by C. Couret). Due to the different altitudes of the two locations, ZSF systematically presents lower O 3 values than ZUG Cooper et al., 2020). To homogenize the two data sets, similar to Cooper et al. (2020), we compared simultaneous observations carried out at ZUG and ZSF from 2002 to 2010. To minimize the possibility that thermal wind circulation and transport from PBL affected our analysis, we calculated the average differences between monthly nighttime (00:00-04:00 UTCþ1) values at ZUG and ZSF ( Figure  S2): In good agreement with Cooper et al. (2020), ZSF observations were lower than ZUG by 1.21 ppb, on average. Then, the O 3 values at ZSF were increased by 1.21 ppb to avoid introducing artificial discontinuities related to the different sampling locations.

Near-surface O 3 in Europe
To obtain information on the variability of near-surface O 3 over specific regions of the European continent, we considered the AirBase (version 7) data set compiled by the European Environmental Agency (2020). AirBase version 7 provides hourly O 3 mean values from 1996 to 2012 for each single air quality station included in the database. We aggregated the measurement sites in five geographical regions (hereinafter "source" regions) which, according to the results reported by Section 2.2, were characterized by high emission sensitivities for CMN, JFJ, and ZUG/ZSF (Figures 1 and 2 Figure S3. The number of reporting stations varied across the years as a function of each "source" region. As an instance, for the NI region, we recorded a minimum of 20 stations during 1996-1998 and a maximum of 200 stations in 2010. For each "source" region, we calculated the time series of the O 3 monthly mean values and monthly percentiles (5 th , 25 th , 50 th , 75 th , and 95 th ). First, for each single measurement station, we calculated the time series of the O 3 monthly mean values. For each "source" region, these values were averaged on a monthly basis to obtain the regional time series of monthly mean values and percentiles. These monthly values were further averaged for obtaining the seasonal (i.e., DJF, MAM, JJA, and SON) averaged values of mean O 3 and percentiles over each single "source" region.

Trend and multi-annual O 3 analysis
To assess the long-term O 3 tendency at CMN and to make a comparison with JFJ and ZUG/ZSF, we applied different methodologies for different data selection. For each measurement site, the trend estimation was performed for the observed monthly mean O 3 as well as for the time series of anomalies which should provide also information about the impact of climate variability on interannual O 3 variations (Logan et al., 2012;Oltmans et al., 2013;Lin et al., 2014Lin et al., , 2015Tarasick et al., 2016;Cooper et al., 2020). The time series of O 3 anomalies were calculated by subtracting the average O 3 monthly values over the full period 1996-2016. A data coverage threshold of 75% was imposed for the calculation of the monthly averages.
Since the focus of our study are O 3 time series at three high-mountain sites, we cannot neglect the possible impact of the air mass transport from the PBL. To take Art. 8(1), page 4 of 27 Cristofanelli et al: Analysis of near-surface ozone in the northern Mediterranean basin into account the impact of different atmospheric mixing conditions to O 3 diurnal cycles, we categorized O 3 data as a function of the time of day. Nighttime (00:00-04:00 UTCþ1) data were selected as less influenced by transport of air masses from PBL, while daytime (10:00-18:00) data were considered more impacted by PBL air masses. The whole data set (00:00-23:00) was also considered for comparison. We adopted narrower time windows in respect than the strategy used in the Tropospheric Ozone Assessment Report (TOAR; e.g., Schultz et al., 2017) to reduce the possibility that the results could be biased by the presence of data representative of the transition from the more to the less PBL-affected time periods.
To estimate the long-term O 3 trend, in agreement with Cooper et al. (2020), we applied a linear regression model to observed monthly mean O 3 and to the time series of anomalies for each measurement site. As reported by Cooper et al. (2020), this method is appropriate for long-term time series, and it is robust against outliers when at least 20 years of data are considered (which is our case, indeed). Moreover, it takes into account autocorrelation. The linear regression model can be summarized as follows: Where y is the monthly mean values of observed O 3 or anomaly, t is the monthly index from January 1996 to December 2016, a is a constant, b is the linear trend, g and d are coefficients describing the 12-month harmonic series of the seasonal cycle (with M ¼ 1, . . . , 12), and R t describes the contributions from the autocorrelation (with length ¼ 1) and a normal random error.
With the aim of integrating the results from the linear regression model, other statistical methods have been used for calculating the long-term trend at CMN, JFJ, and ZUG/ZSF. The Theil-Sen estimator (Theil, 1950;Sen, 1968) is known to be a suitable trend metric even with nonnormal data to be resistant to outliers (Carslaw and Ropkins, 2012;Lefohn et al., 2018). Given a set of n(x, y) pairs, the slopes between all pairs of points are calculated: The Theil-Sen estimate of the trend slope is the median of all these slopes. This metric was used in the first TOAR (see Lefohn et al., 2018) for estimating the magnitude of O 3 trends. One assumption of the usage of the Theil-Sen estimator is that the time series should be monotonic. This is clearly not the case for CMN, which shows an O 3 increase in 2004-2008. However, Theil-Sen estimator can be safely used for JFJ and ZUG, allowing a comparison with results from the linear regression model. For the Theil-Sen slope, a 95% confidence interval was also provided using a bootstrap simulation (Carslaw and Ropkins, 2012). Finally, we consider the Mann-Kendall estimator for calculating the P-value of trends. Similar to Theil-Sen, the Mann-Kendall test was used in the framework of TOAR : It is also resistant to outliers; it does not require assumptions regarding functional form or statistical distribution for the data.
The investigation of the long-term trends using the three mentioned approaches will allow to better insert our results in the perspectives provided by international reference initiatives for the assessment of tropospheric O 3 trends like TOAR and the recent paper by Cooper et al. (2020).

Multi-annual O 3 variability by seasonal and trend decomposition using loess (STL)
To better compare the multi-annual variability of O 3 at the three measurement sites, we analyzed the time series using the "STL" methodology, which decompose a time series into seasonal, trend, and irregular components (Cleveland et al., 1990). Since the STL results are particularly sensitive to the definition of the smoothing parameter to fit the seasonal component, the STL was run in two different configurations: (1) by replacing the seasonal smoothing by taking the mean seasonal cycle and (2) by setting the seasonal trend smoothing parameter to 7 (as recommended by Cleveland et al., 1990). Moreover, two further runs have been carried out by adopting (or not) a "robust" fitting in the Loess procedure (Figures S4-S9). As indicated by Cleveland et al. (1990), the STL "robust" estimation is needed when data checking indicates non-Gaussian behavior in the time series (i.e., presence of "outliers"). This can be the case for CMN since the July 2006 monthly value outstands like an "outlier" from the population of the monthly mean values. The use of the STL methodology does not provide a "direct" quantification of long-term trend but allows to compare the multi-annual O 3 variation at the three measurement sites, thus helping in better contextualize CMN observations by comparison with other high-altitude "reference" observations in Europe.

Air mass transport analysis
Air mass transport analysis to the receptor sites is based on backward simulations with Lagrangian particle dispersion model FLEXPART (Stohl et al., 1998(Stohl et al., , 2005, driven by operational 3-hourly meteorological data at 1 Â 1 resolution from the European Centre for Medium-Range Weather Forecasts. FLEXPART model calculates the trajectories of tracer particles using the mean winds from the analysis fields, random motions describing turbulence (Stohl and Thomson, 1999), and convection. From each measurement station, 40,000 particles were released every 3 h from 1996 to 2016 and followed for 20 days backward in time. The aim of our model analysis is to investigate the transport pattern influencing the monitoring sites from the regional to the intercontinental scales.
Using 20 days backward trajectories should be long enough to capture the transport from the most relevant source regions. Moreover, the choice of simulation length was motivated by the fact that the enhancement value of every additional simulation day decreases rapidly with time backward. Indeed, due to turbulent mixing and convection, the emission contributions from various regions become more and more well mixed and start forming the baseline, furthermore the model errors growth with time.
Following the consolidated approach by Seibert and Frank (2004), FLEXPART calculations allow to obtain the sensitivity of the receptor sites to geographical sources, that is, the source receptor relationship (SRR). For a particular grid cell in the spatial domain, the SRR is proportional to the particle residence time in that cell and measures the simulated mixing ratio that a source of unit strength (1 kg s -1 ) in the cell would produce at the receptor sites (Stohl et al., 2009). The SRR is hence a proxy for the influence that chemical species emissions from individual regions may have on the observation at the site. To investigate the potential impact of PBL air masses, we considered the SRR related to the vertical layers (0-100 m a.s.l.). To take into consideration air masses originating at the interface between PBL and the free troposphere, the layer 100-3,000 m a.s.l. was investigated. In general, no significant differences have been found for the two layers (here not shown). For this reason, in the following, only the results for the lowermost atmospheric layer are showed: In case significant deviations occurred between the two layers, these will be specifically discussed.
Besides absolute SRR, we also calculated SRR anomalies by comparing SRRs obtained for specific seasons with those calculated for a reference period. We used these analyses to assess the different sensitivities of CMN, JFJ, and ZUG/ZSF to European source regions as well as the existence of anomalies in the transport regime at in connection with the occurrence of detected O 3 anomalies.
To study the different sensitivity with respect to the surface emissions, for the "reference" period 1996-2016, we calculated the seasonally averaged SRRs for the "surface" layer 0-100 m a.s.l. (Figure 1). The resulting patterns vary among the different seasons, and they are in agreement with previous analysis by Maione et al. (2008), revealing that CMN generally appears to be more sensitive to emissions occurring over the European domain than JFJ. Figure 2 reports the seasonal averaged differences of SRR between CMN and the other sites expressed in % (positive values denote higher sensitivity for CMN). It is evident that, with respect to JFJ and ZUG/ZSF, CMN is more sensitive to potential emissions occurring over northern Italy, Mediterranean Sea, northern Africa, and continental Europe (central and eastern Europe, especially). JFJ and ZUG/ZSF appear to be more sensitive to potential emissions occurring over the European NW quarter (e.g., France).

Long-term near-surface O 3 variability
Long-term O 3 time series at CMN, JFJ, and ZUG/ZSF are presented in Figure 3 as monthly mean values calculated based on the 1-h averages with a 75% data coverage criterion. Absolute mixing ratios are in the same range for all the stations (overall means: 53.3 + 1.1 ppb at CMN; 52.5 + 0.8 ppb at JFJ, and 50.4 ppb + 0.8 ppb at ZUG/ZSF; + 95% confidence level), and the seasonal cycles are in Art. 8(1), page 6 of 27 Cristofanelli et al: Analysis of near-surface ozone in the northern Mediterranean basin phase with spring-summer maxima and winter minima, but the mean annual peak to peak amplitude is larger at CMN (average value: 21.1 ppb) than at JFJ and ZUG/ZSF (14.9 and 14.7 ppb). The lower O 3 values observed at ZUG/ ZSF with respect to JFJ can be explained considering the lower altitude of this measurement site . Conversely, the higher O 3 values observed at CMN with respect to ZUG/ZSF can mimic the larger exposure of this Mediterranean site to photochemical O 3 production and air mass transport from PBL. This is also testified by the larger deviations of O 3 at CMN during warm months with respect to JFJ and ZUG/ZSF, suggesting that during spring-summer CMN is more affected by photochemically produced O 3 . All-time series were also deseasonalized by using a robust Loess smoothing by setting the seasonal trend smoothing parameter to 7 (as recommended by Cleveland et al., 1990). The time series of deseasonalized monthly averages (Figure 3, middle panel) reveal persistent patterns at both stations (like the smooth harmonic-like pattern between 1996 and 2003) and similar O 3 values in the most recent years (since 2009). The Pearson correlation coefficient between deseasonalized monthly mean values for CMN and both the Alpine stations is 0.59 + 0.09 (95% confidence level). Table 1 reports the Pearson correlation coefficient for each single season and for different data selections (all data, nighttime, daytime): Positive correlations characterized CMN and ZUG/ZSF during all the seasons but not in autumn for JFJ. In general, the Pearson coefficient is lower for the daytime data selection, probably indicating the larger influence of PBL air masses at the lower altitude CMN site.
The average differences in the monthly mean O 3 values for JFJ and ZUG/ZSF with CMN (positive DO 3 indicated higher values at CMN with respect to the reference Alpine stations) for the whole period were 0.7 + 0.5 ppb (P ¼ 0.05) and 3.1 + 0.5 ppb (P ¼ 0.05), respectively. CMN reports consistently higher values with respect to JFJ and ZUG/ZSF between 2004 and summer 2008, with DO 3 averaging 4.3 + 1.0 ppb (P ¼ 0.05) and 6.5 + 0.6 ppb (P ¼ 0.05), respectively. When looking into the whole 1996-2016 data set, other periods showed comparable differences between CMN and JFJ or ZUG/ZSF (e.g., during 1997, 2011, or 2012) but with a more limited extension in time.

Long-term O 3 tendencies/trends
In the attempt of describing the long-term O 3 tendencies (and trends) at CMN in the context of others European high-mountain stations, here we reported (Table 2) the results for the different trend metrics (linear model, Theil-Sen, and Mann-Kendall) for the three measurement sites as well as for the different data selections (all data, daytime, and nighttime). Moreover, analyses are reported for the O 3 monthly anomalies as well as for the original O 3 monthly means (Table S1). The magnitude of the tendencies/trends (as provided by the linear fitting model and the Theil-Sen estimator) and their confidence interval/P values (as provided by the linear fitting model, Mann-Kendall, and bootstrap analysis for Theil-Sen) appear to be rather consistent among the data selections (all data, nighttime, or daytime) and the inspected variable (anomalies vs. actual observations). All the measurement sites reported negative tendencies over the period 1996-2016, with the highest magnitude being observed at CMN and the lowest at JFJ. The slope values reported for CMN appeared to be about 2 times higher than those calculated for ZUG/ ZSF and 3-4 times higher than for JFJ. Even considering the slightly different observation periods  as well as the different data selection (from 20:00 to 7:59 UTCþ1), our results fit nicely with the recent work by Cooper et al. (2020) for ZUG/ZSF (-0.08 ppb yr -1 with P ¼ 0.00), with O 3 that appeared relatively unchanged at JFJ. However, despite Cooper et al. (2020), who reported a weak positive tendency from 1995 to 2018 (þ 0.02 ppb  (Table 3): May-September (when thermal transport and PBL mixing is higher) and November-February (when thermal transport and PBL mixing is minimized). For all the measurement sites, negative trends are stronger during the warm months. In particular, the daytime data selection reports an evident negative trend from May to September (MJJAS) also for JFJ. The robust negative trends during the warm months are consistent with the results by Cooper et al. (2020), while Gaudel et al. (2018) reported no robust negative tendencies for JFJ and ZUG/ ZSF over 2000-2015 for JJA. In our analysis, the stronger negative trends during the daytime of warm months would point toward a decrease in the O 3 within air mass from the European PBL. Based on these results, the annual negative trends at CMN and ZUG/ZSF (as well as the weaker negative tendency at JFJ) reported in Table 3 appeared to be mostly driven by the O 3 decrease observed during MJJAS. Despite the other measurement sites, CMN showed robust negative trends even during the cold months.
With the purpose of investigating the reasons behind the stronger negative trends observed at CMN, we used the STL decomposition to visualize the multi-annual O 3 fluctuations at the three different sites. Figure 4 reports the results of STL robust decomposition for the original O 3 monthly data sets at CMN, ZUG/ZSF, and JFJ using the smooth parameter t ¼ 7. Results from other STL runs can be found in the supplementary material. In general, the "trend" components appeared to be consistent among the different STL runs for each measurement site, but some differences can be pointed out. As an instance for CMN, a lower peak around 2008 was obtained if compared with the STL using the "not robust" approach and t ¼ "seasonal," while for JFJ and ZUG/ZSF, a lower peak in 2003 can be observed. We are particularly interested in the "trend" component as it depicts the multi-annual variability of the time series. In agreement with the trend analysis, all the time series showed generally lower values of the "trend" component in the most recent period: The average mean value of the last 5 years (i.e., January 2011 to December 2016) is 51.0 + 0.2 ppb at CMN, 51.7 + 0.1 ppb at JFJ, and 49.7 + 0.1 ppb at ZUG/ZSF. In the previous 15 years, the 5-year average mean value ranged from 54.5 + 1.3 ppb to 53.6 + 1.9 ppb at CMN, from 52.9 +   Figure 3). Hereinafter, these periods will be referred to as "persistent" O 3 episodes.
With the aim to better characterize the "persistent" O 3 episodes, for each season of the year, we calculated the seasonal probability density functions (PDFs) of nearsurface O 3 at CMN, JFJ, and ZUG/ZSF for the whole data set and for the selected episodes ( Figure 5). PDFs are calculated on daily (i.e., average mean values from 00:00 to 23:00 UTCþ1) O 3 data.
In DJF (i.e., December-January-February), the persistent O 3 episodes were related to a shift of the PDF peak at CMN (by about 9 ppb) toward higher values with a left-hand We argued whether the presence of these "persistent" O 3 episodes concurred to explain the stronger decreasing O 3 trends observed at CMN with respect to JFJ and ZUG/ZSF. To make a rough test about the impact of the "persistent" O 3 episodes to the CMN trend, we substituted the actual O 3 observations from March 2004 to May 2008 with an artificial time series taking into account the typical O 3 seasonal cycle (calculated with "periodic" STL decomposition, see Figure S10). By applying the linear fitting model to the O 3 anomalies, we obtained a trend of -0.18 + 0.05 ppb  (P < 0.01). The obtained trend value is still comparable to the value obtained for the original time series; thus, we can argue that the "persistent" O 3 episodes were not the principal reasons of the stronger O 3 trend at CMN. Even if the "persistent" O 3 episodes in 2004-2008 appeared as not strongly influencing the long-term trend calculation at CMN, it is important to understand the underlying processes that determined them.
Several reasons can be responsible to explain the O 3 features at CMN in 2004-2008: anomalies in air mass transport regimes on different spatial scales (from global to regional), variability of precursor emissions, or in meteorological conditions affecting O 3 production/removal (e.g., prolonged heat waves or rainfall anomalies). In the following, we provided a detailed analysis of a subset of possible processes that can impact O 3 variability at CMN. In particular, we compared the diurnal O 3 variability at CMN and at the rural site FEB, located in the Apennines foothills, to assess the possibility that local thermal transport or PBL mixing occurring at diurnal scales could   No evident change in the diurnal cycle is observed for DJF between persistent O 3 episodes and normal periods. However, the average seasonal diurnal cycles show an overall shift toward higher values during the positive anomaly years for MAM, JJA, and SON. This appears to be consistent with the CMN observations. Higher O 3 values were observed at FEB during the central part of the day, when upward air mass transport may reach mountain regions and photochemistry is more active (Cristofanelli et al., 2015, and reference therein). With respect to normal periods, the amplitudes of the average diurnal cycles increased in JJA 2005JJA -2007 and SON 2006and SON -2007, thus indicating a possible enhanced impact of processes occurring at diurnal scale (i.e., thermal air mass transport and photochemistry): Local effects coherently appear to be more relevant at FEB than at CMN, likely due to the relative distance of the two sites from the closer pollutant source region (i.e., the Po basin).

Air mass transport analysis
FLEXPART-based SRRs (Section 2.5) were analyzed to assess a possible influence of the air mass transport regime on the observed O 3 episodes. Due to the lower systematic exposure to direct transport of PBL air masses (compared to ZUG/ZSF), here we only considered JFJ as a reference site to investigate the anomalies observed at CMN.
With the aim of exploring the possible relationship between O 3 variability and air mass transport at CMN and JFJ, we calculated the SRR tagged with specific geographical "source" regions (i.e., SRR reg ). The SRR reg was calculated for spatial domains that are well known to be characterized by high O 3 (see, e.g., Gaudel et al., 2018): "northern Italy" (including the Po basin, "NI"), "western Europe" (WE), "continental Western and Eastern Europe" (CWE and CEE, respectively), and "Mediterranean basin" (MED). Moreover, we also set three geographical boundaries to investigate the possible contribution of longrange transport: "North Atlantic" (NA), "North America" (NAM), and the rest of the Northern Hemisphere (NH), see Figure 8. On average, these source regions accounted for 97.6% of the global SRR at CMN and 92.6% at JFJ for the layer 0-100 m a.s.l. (the same values were found for the layer 100-3,000 m a.s.l.). Figure 9A and B provides indications about the relative importance of each source region in determining the total SRR at CMN and JFJ for the different seasons. Due to their wide catchment areas, the "long-range" transport regimes (i.e., NH, NA, NAM) accounted for 56.2% (CMN) and 66.7% (JFJ) of total SRR in all the seasons, with a prominent contribution of NH and NA. A seasonal variation occurred for the relative importance of source regions. For CMN, NI and CWE are maximized in summer (11% and 17%, respectively). As expected, MED contribution appeared to play an important role in affecting atmospheric variability at CMN with fractional contribution ranging from 9% in winter to 11% in autumn. The contributions by the "long range" transport regime NH appeared to be strongly minimized in summer. For JFJ, the contributions related to WE and CWE regions are maximized in summer (11% and 13%). As for CMN, also at JFJ, the MED contribution is maximized in autumn (8%), while the contributions from "long range" NH is minimized in summer. It is possible that the contribution from NI (for CMN) and WE (for JFJ) can be underestimated due to the uncertainty in the simulation of vertical transport in this region characterized by complex orography (Zhang et al., 2017). Figures 10-13 report the time series of seasonal SRR reg /SRR tot for the different source regions at CMN (red) and JFJ (blue). Here, the ratio SRR reg /SRR tot is used as a proxy of the relative importance of each single source region in potentially affecting atmospheric observations at the measurement sites. For both the measurement sites, the seasonal SRR reg /SRR tot values show evident interannual fluctuations. However, decreasing tendencies in SRR reg /SRR tot can be observed for the source regions related to "long-range" regimes at JFJ. As supported by the application of the nonparametric Mann-Kendall test (Mann, 1945), this is especially evident for the source regions NA in spring (-0.8% yr -1 , with P < 0.05) and autumn (-0.9% yr -1 , P < 0.01), NAM in summer (-1.4% yr -1 , P ¼ 0.05) and autumn (-1.8% yr -1 , p < 0.01), NH in winter (-0.6% yr -1 , p < 0.10) and autumn (-0.8% yr -1 , P ¼ 0.10). On the contrary, an increasing tendency can be observed for the contributions related to "regional" source regions: NI in spring (þ1.9% yr -1 , P ¼ 0.10) and autumn (þ4.2% yr -1 , P < 0.01); WE in winter (þ2.1% yr -1 , P < 0.05), spring (þ1.1% yr -1 , P < 0.10), and autumn (þ2.3% yr -1 , P < 0.01); CWE in spring (þ1.2% yr -1 , P < 0.05) and autumn (þ1.5% yr -1 , P < 0.05); and CEE in spring (þ5.7% yr -1 , P < 0.05) and autumn (þ5.0% yr -1 , P ¼ 0.11). At CMN, long-term tendencies of SRR reg /SRR tot are less evident: An increasing tendency can be observed for the source regions MED (in autumn, þ2.3% yr -1 P < 0.10) and NAM (in winter, þ1.4% yr -1 , P < 0.05).
To identify possible drivers of the increase of O 3 at CMN with respect to JFJ over the period 1996-2016, we calculated the linear correlation matrix between seasonal DO 3 , O 3 at CMN, O 3 at JFJ, and the SRR reg /SRR tot values for the different source regions at both the measurement sites ( Figures S11-S14). The first outcome of this analysis is that for all the seasons, the interannual variability of DO 3 was mainly related to the variability of O 3 at CMN (ranging from R ¼ .76, in JJA to R ¼ .96 in SON), that is, the variability in the deviations of O 3 between CMN and JFJ were mostly due to the O 3 variability at CMN.
Less clear results were obtained for the relationship between DO 3 and SRR reg /SRR tot values. In winter, the strongest relationships were found between DO 3 and SRR reg /SRR tot related to the source region CWE at JFJ (R ¼ -.051) and the source region NAM at JFJ (R ¼ 0.36). For spring, a positive correlation (R ¼ 0.69) was found with SRR reg /SRR tot related to the source region NI at CMN, while negative correlation was found with source region MED at CMN (R ¼ -0.48). For summer, a positive correlation was found for the source regions NH (R ¼ 0.65) and MED (R ¼ 0.36) at CMN and JFJ (R ¼ 0.61), while a negative correlation was found with the source region NI at JFJ (R ¼ -0.48) and CWE at CMN (R ¼ -0.68) and JFJ (R ¼ -0.51). In autumn, a positive correlation was tagged with SRR reg /SRR tot related to the source regions NI at CMN (R ¼ 0.51) and NAM at JFJ (R ¼ 0.45), while a negative correlation was found with the source region CWE at JFJ (R ¼ -0.62). Overall, the positive correlations between DO 3 and SRR reg /SRR tot related to European source regions at CMN (i.e., NI in MAM and autumn; MED in summer) can be explained in terms of advection of air masses enriched in photochemically produced O 3 from the PBL. The negative correlation related to the source regions MED in spring and CWE in summer can be tentatively explained by the occurrence of other drivers with respect to air mass transport (e.g., impact of mineral dust or meteorological variability on tropospheric O 3 , see Pausata et al., 2012;Duchi et al., 2016). On the other hand, because DO 3 variability is mostly driven by O 3 at CMN, a clear explanation for the relationship between DO 3 and SRR reg /SRR tot at JFJ (i.e., CWE and NAM in winter; MED, NI, and CWE in summer; and NAM and CWE in autumn) cannot be given.
At both the measurement sites, the SRR reg /SRR tot time series are characterized by a strong variability. For specific source regions, peaks or minima in the SRR reg / SRR tot were evident during years and seasons for which the persistent O 3 episodes were observed at CMN. In DJF 2007, lower than average SRR reg /SRR tot values characterized CMN and JFJ for the source region CWE. Moreover, in DJF 2008, low SRR reg /SRR tot can be observed at CMN for the source region CEE ( Figure 10). Thus, the high O 3 values at CMN can be possibly related to a decreased contribution of continental PBL air masses depleted by NO titration. However, it should be noted that similar low SRR reg /SRR tot values were observed for other years (i.e., 2012(i.e., and 2014(i.e., for CWE and 2004(i.e., , 2014(i.e., , 2016  than average SRR reg /SRR tot were tagged to the MED region. This is consistent with a possible increase of near-surface O 3 due to a lower occurrence of mineral dust transport from northern Africa, as also supported by Duchi et al. (2016) and Bauer et al. (2004). However, a similar anomaly in SRR reg /SRR tot was evident for MED also at JFJ. For JJA (Figure 12), high SRR reg /SRR tot affected the CMN source regions CEE (in 2006), NA (in 2004), and NH (2005-2007: It can be argued that variability in regional (CEE) and long-range (NH and NA) transport can partially concur to the observed anomaly. These anomalies in air mass transport were less evident at JFJ.
For SON 2006 (Figure 13), high SRR reg /SRR tot values characterized the source regions NI, NA, and NAM at CMN, while the highest SRR reg /SRR tot value is observed for NH in SON 2007. Similar anomalies were also observed for JFJ; it is thus difficult to definitively attribute the observed O 3 anomalies at CMN to these anomalies in the air mass transport.
Further hints for a possible attribution of the "persistent" O 3 episodes at CMN can be obtained by the analysis of the spatial distribution of seasonal SRR over the European continent ( Figure 14). Negative (positive) anomalies in the SRR values are evident for CMN (JFJ) over NI during DJF 2007 suggesting that changes in the regional scale atmospheric circulation could affect O 3

Near-surface O 3 analysis over Europe
To explore the possibility that "anomaly" in near-surface O 3 over the European continent could contribute to the "persistent" O 3 episodes at CMN, we considered the multiyear fluctuations of near-surface O 3 over the "source" regions in the European domain: WE, CWE, CEE, NI, and MED. To this aim, for each "source" region, we calculated the seasonal O 3 anomalies (i.e., the differences between actual seasonal mean values and seasonal cycle averaged over the whole observation periods) obtained from the monthly mean values of the Airbase data set. We decided to use this metric to magnify the interannual variability of O 3 over the "target" regions in a very effective way. To obtain specific information about the occurrence of elevated O 3 episodes, the anomalies for 5 th , 25 th , 50 th , 75 th , and 95 th percentiles were also calculated. Figure  S15 in the supplementary material reports the time series of monthly O 3 mean values and percentiles for each European "source" region, while Figures S16-S20 report the time series of seasonal anomalies. In the following, we analyzed the seasonal difference of O 3 between CMN and JFJ (again considered as a "reference" station) as a function of seasonal SRR reg /SRR tot and seasonal O 3 anomalies over each "source" region ( Figure 15). This analysis is devoted to provide, for each "source" region, an overview of the possible impact of transport and near-surface O 3 anomalies to the occurrence of enhanced (with respect to JFJ) seasonal O 3 values at CMN. Due to the limited temporal For winter (DJF), the highest DO 3 values were related to negative anomalies of near-surface O 3 over all the European "source" regions. As already pointed out in Section 3.3.3, the air mass transport contribution from NI and other European continental regions (CWE and CEE) showed among the lowest values for these events. Thus, it can be argued that an anomaly in the air mass transport from the regional PBL can contribute to explaining the O 3 behavior at CMN. During spring (MAM), no clear relation between DO 3 and European near-surface O 3 was observed for the "persistent" O 3 episodes, with positive and negative anomalies of the European near-surface O 3 equally distributed among the persistent O 3 episodes. The same is true for summer (JJA). However, in agreement with the results in Section 3.3.3, a tendency for observing high DO 3 with high (low) SRR reg /SRR tot over NI (MED) during spring is evident. As deduced by the high DO 3 , enhanced seasonal O 3 at CMN was related to high near-surface O 3 over all the European "source" regions in autumn (SON). For NI, these positive nearsurface O 3 anomalies were also related to high air mass contribution (as denoted by the high SRR reg /SRR tot values). Similar results are obtained when the seasonal 95 th percentiles of near-surface O 3 over the European "source" regions are considered (see Figure S21 in the supplementary material).

Discussion and conclusions
In this article, we investigated the long-term variability of near-surface O 3 at CMN, a high-altitude site (2,165 m a.s.l.) in the Mediterranean region (Italy, northern Apennines), where continuous measurements of this short-lived climate forcer exist since 1996. In particular, we investigated the multi-annual variability in comparison with the highly mature data sets from two Alpine measurement sites: Jungrfaujoch (Swiss Alps) and Zugspitze (German Alps Based on the analysis of average diurnal O 3 variability at CMN and FEB, processes occurring at diurnal scale (i.e., thermal transport of air masses or local photochemistry) appeared to have played only a minor role in the occurrence of the positive O 3 anomaly at CMN. Even if not conclusive, our analyses suggested that variability in the regional-scale transport can partially explain the appearance of O 3 anomalies at CMN. During the period 1996-2016, we found correlations between DO 3 (i.e., the difference between O 3 at CMN and at JFJ) with transport of air masses from both continental and long-range source regions to CMN and JFJ. However, for the "persistent" episodes in 2004-2008, we were able to detect significant variations of SRR reg /SRR tot only for specific seasons, years, and source regions. Even if we are not able to provide Figure 15. Analysis of the relationship between the seasonal difference of O 3 between CMN and JFJ together with air mass transport contributions from the European "source" regions and the seasonal near-surface O 3 anomalies over these regions. A further point of interest is related to the long-term variability of SRR reg /SRR tot values for some selected regions. In particular, for JFJ, we pointed out a decreased influence related to long-range transport during the more recent observations (2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016) compared to the early years (1996)(1997)(1998)(1999)(2000). This was accompanied by a relative increase of SRR reg /SRR tot values from regions in the European domain. Similar coherent long-term SRR reg /SRR tot variability cannot be detected for CMN.
All these points only partially reconcile the existing deviations between historical time series of near-surface O 3 at CMN and JFJ in terms of variability of atmospheric transport (and different sensitivities of the measurement sites to specific "source" regions). With the purpose of evaluating the role of emission of O 3 precursors or the role of processes occurring in the PBL (e.g., NO titration), only the "surface" layer (0-100 m a.s.l.) was discussed in this work. The limitation of the model in reproducing mesoscale circulation in complex mountain terrain can be a caveat of our analysis, especially for the summer periods. It is thus possible that contributions related to "special" regional episodes like heat waves can be underestimated.
Although the underlying processes driving the detected changes in the regional-scale circulation diagnosed by FLEXPART need to be identified, these results emphasize that several factors, with a special emphasis on atmospheric transport regime and variability of O 3 over specific "source" regions, must be considered when long-term O 3 variability is discussed. Followup studies may be related to a detailed analysis of longterm tendencies of SRR reg /SRR tot .

Data accessibility statement
The near-surface O 3 data series (1-h average values) used in this work for CMN (1996), JFJ (1996 and ZUG/ZSF (1996 are available in the archive the World Data Centre for Reactive Gases (http://ebas.nilu.no). ZSF data for year 2011 are directly available by C. Carnout. Furthermore, CMN data can be accessed by the MOVIDA-Multistation web system (http:// shiny.bo.isac.cnr.it:3838/plot-multistats-en/). The nearsurface O 3 data series for FEB are available by ARPAE Emilia-Romagna by their web site https://www.arpae. it/v2_rete_di_monitoraggio.asp?p¼RE&idlivello¼1637 (last accessed July 2020). The AirBase (version 7) data set compiled by the European Environmental Agency is available at https://www.eea.europa.eu/data-and-maps/data/ airbase-the-european-air-quality-database-7 (last accessed July 2020).

Supplementary files
The supplemental files for this article can be found as follows: Figures S1-S21. Table S1. Docx