**E**ff**ects of Ozone and Clouds on Temporal Variability of Surface UV Radiation and UV Resources over Northern Eurasia Derived from Measurements and Modeling**

**Natalia E. Chubarova 1,\* , Anna S. Pastukhova 1, Ekaterina Y. Zhdanova 1, Elena V. Volpert 1, Sergey P. Smyshlyaev <sup>2</sup> and Vener Y. Galin <sup>3</sup>**


Received: 25 November 2019; Accepted: 23 December 2019; Published: 2 January 2020

**Abstract:** Temporal variability in erythemal radiation over Northern Eurasia (40◦–80◦ N, 10◦ W–180◦ E) due to total ozone column (X) and cloudiness was assessed by using retrievals from ERA-Interim reanalysis, TOMS/OMI satellite measurements, and INM-RSHU chemistry–climate model (CCM) for the 1979–2015 period. For clear-sky conditions during spring and summer, consistent trends in erythemal daily doses (*Eery*) up to +3%/decade, attributed to decreases in X, were calculated from the three datasets. Model experiments suggest that anthropogenic emissions of ozone-depleting substances were the largest contributor to *Eery* trends, while volcanic aerosol and changes in sea surface temperature also played an important role. For all-sky conditions, *Eery* trends, calculated from the ERA-Interim and TOMS/OMI data over the territory of Eastern Europe, Siberia and Northeastern Asia, were significantly larger (up to +5–8%/decade) due to a combination of decrease in ozone and cloudiness. In contrast, all-sky maximum trends in *Eery*, calculated from the CCM results, were only +3–4%/decade. While *Eery* trends for Northern Eurasia were generally positive, negative trends were observed in July over central Arctic regions due to an increase in cloudiness. Finally, changes in the ultraviolet (UV) resources (characteristics of UV radiation for beneficial (vitamin D production) or adverse (sunburn) effects on human health) were assessed. When defining a "UV optimum" condition with the best balance in *Eery* for human health, the observed increases in *Eery* led to a noticeable reduction of the area with UV optimum for skin types 1 and 2, especially in April. In contrast, in central Arctic regions, decreases in *Eery* in July resulted in a change from "UV excess" to "UV optimum" conditions for skin types 2 and 3.

**Keywords:** total ozone content; cloudiness; erythemal radiation; trend; chemical–climate model; ERA-Interim reanalysis; Northern Eurasia; UV resources

#### **1. Introduction**

Ultraviolet (UV) radiation is well-known for its significant influence on human health and the environment. High UV doses have negative effects on skin (erythema (sunburn), skin cancer), and cause eye diseases and immune suppression [1,2]. However, moderate UV doses have positive effects causing vitamin D production [3,4]. Apart from the solar elevation, ozone and cloudiness are the main factors affecting UV level [2,5] and providing significant year-to-year variability of UV radiation. Erythemal

radiation, which is used for characterizing UV health effects, is the UV radiation, weighted by the erythemal action spectrum with maximum efficiency at 280–300 nm and integrated over 280–400 nm [1].

As stated in numerous publications and assessments [5–8], the decrease in stratospheric ozone resulting from anthropogenic emissions of ozone depleting substances (ODSs) into the atmosphere has led to the increase in erythemal radiation at the Earth's surface in the 1980s and 1990s. Since the end of 1990s there is a small increase in total ozone content from 0.3% to 1.2% per decade due to the reduction of the ODSs in the stratosphere but this increase is not statistically significant [8]. However, there is no clear relationship between total ozone column (X) and erythemal radiation for many regions due to the influence of other factors like cloudiness, aerosols, minor gas species, and surface reflectivity [9,10]. The most important atmospheric factor is cloudiness, which significantly affects erythemal radiation and can enhance or mask its trends due to ozone changes.

There are a lot of studies devoted to the temporal and spatial variability of erythemal radiation over the globe or large spatial areas using the retrievals from satellite instruments, and model simulations [9–11], and references therein]. According to these data, different kinds of UV climatologies at various spatial scale have been obtained [12–15]. In addition to the UV climatologies, many publications are devoted to the analysis of temporal changes in *Eery* due to ozone and cloudiness. In [16] the analysis of satellite data over the globe has revealed an increase in erythemal radiation of about 8% at 50◦ S and of about 5% at 50◦ N over the 1978–2008 period due to decrease in X, while the combined effect of ozone and cloudiness led to the increase of about 2% at 50◦ S and 6% at 50◦ N. The changes in erythemal radiation, evaluated over the 1979–2010 period using TOMS/OMI retrievals, have also shown positive trends up to 5% at 45◦–50◦ S in spring–summer seasons [17]. A statistically significant decrease in Lambertian Equivalent Reflectivity due to clouds (3.6 ± 0.02%) over the 1979–2011 period was derived from different satellite measurements [18,19]. However, satellite data have a large uncertainty due to problems with accounting for absorbing aerosol and due to the restrictions in UV retrievals over bright snow/ice surfaces [20,21] compared with direct ground-based measurements.

Many publications are devoted to the analysis of temporal variations in UV radiation obtained over local sites ([2] and references therein, [22–33]). According to these publications, an increase of UV radiation has been revealed over many European sites since the middle of 20 century, mainly due to ozone and cloud effects [28,29,31,33]. Over several locations (e.g., Thessaloniki) the increase in *Eery* to some extent can be also attributed to the decrease in aerosol optical thickness [34]. However, in the west Arctic area at Hornsund (77◦00 N, 15◦33 E) no pronounced trend in *Eery* has been found since 1983 due to complicated changes in cloud cover, which is one of the basic drivers of the long-term UV changes over this region [30].

Over Europe the UV climatology and long-term variability in erythemal radiation was analyzed within the framework of the COST Action 726 project [32] using accurate model simulations, the results of ERA-40 reanalysis over the five-decade period since 1950, and ground-based measurements. These data also demonstrate pronounced variations in *Eery* due to ozone and cloudiness with some specific features over different areas.

For evaluating variations in *Eery* in the past and future, global Earth-system models from the CMIP-5 (Coupled Model Intercomparison Project Phase 5) experiment, the ensemble mean of the Community Earth System Model version 1 (CESM1) and the Whole Atmosphere Community Climate Model (WACCM) with interactive chemistry ensemble of model ozone predictions have been used for estimating the UV variability over the 1955–2100 period [5]. Using this approach, the effects of various factors on future UV variations were estimated. Chemistry–climate models (CCMs) have also been used for estimating changes in UV radiation in the past and future [9,11,35–38]. In [35] the projections of erythemal irradiance from 1960 to 2100 have been made using radiative transfer calculations and projections of ozone, temperature and cloud change from 14 CCM, as part of the CCMVal-2 activity of SPARC (Stratosphere-troposphere Processes And their Role in Climate) project. According to this study annually mean erythemal radiation in the 2090s relative to 1980 will be on average approximately 12% lower at high latitudes in both hemispheres, and 3% lower at midlatitudes. At northern high latitudes, the increase in cloudiness towards the end of the 21st century will reduce the erythemal irradiance by 5% with respect to the 1960s.

In Ref. [37] the analysis of *Eery* variations has been made by using the clear-sky data from the first phase of the Chemistry–Climate Model Initiative (CCMI) as input to the TUV (Tropospheric Ultraviolet and Visible) radiation model. This analysis projects an increase in average *Eery* of about 2–4% in 2100 in the tropical belt (30◦ N–30◦ S) and a 1.8% to 3.4% increase in the midlatitudes in the Southern Hemisphere for RCP (Representative Concentration Pathway) 2.6, 4.5 and 6.0, compared to 1960s, which partly contradict to the results obtained in [35]. The projected increase in *Eery* reported in [37] results from the assumption that the atmospheric aerosol loading will decrease greatly over the course of the 21st century, which is debatable. The analysis of erythemal radiation according to the CCMI simulations [2] projects *Eery* to decrease by 5–15% in the northern hemisphere during summer and autumn mainly due to ozone recovery in 2085–2095 compared with 2010–2020 according to the RCP 6.0 scenario. However, some other factors like cloudiness, aerosol and surface reflectivity can also play a significant role over several regions. In [39] negative *Eery* changes over Northern midlatitudes of about 6–8% due to ozone and clouds were found by the end of 21 century, which are in agreement with [35].

For a wide range of biological and medical studies, it is important to interpret long-term changes in UV radiation from the point of health effects. The approach proposed in [14] distinguishes health effects between harmful and beneficial using the term "UV resources". The UV resources characterize UV radiation from the point of beneficial (vitamin D production) or adverse (sunburn) effects on human health. Their classification is given in the categories of UV deficiency, UV optimum and UV excess and takes into account the skin type and the fraction of the skin surface that is exposed to UV radiation. A similar approach has been previously developed in [40]. Changes in UV resources over the 1979–2015 period were first simulated for Moscow conditions in [28], where a transition from the UV optimum to UV moderate excess conditions for the population with the most vulnerable skin type I has been obtained in spring.

The main aim of this study is to identify temporal variability in Eery and the change in UV resources due to ozone and cloudiness between 1979 and 2015 based on the TOMS/OMI satellite datasets, and UV retrievals from the reanalysis and a chemistry–climate model over Northern Eurasia (40◦ N–80◦ N, 10◦ W–180◦ E). We chose this period since it was characterized by the most pronounced changes in *Eery*, and the quality of data is the most reliable. In the analysis we used the updated TOMS and OMI satellite UV retrievals using the new Macv2 absorbing aerosol correction [41], as well as the ERA-Interim Reanalysis data [42]. For revealing the causes of changes in ozone we applied model runs of the Russian Chemistry–Climate Model (CCM), which was jointly developed by the Institute of Numerical Mathematics RAS (INM) and Russian State Hydrometeorological University (RSHU) [43,44].

#### **2. Methods and the Data Description**

Long-term variability in *Eery* over Northern Eurasia was estimated using different datasets, which include the UV data from model, reanalysis, and satellite measurements. For evaluating the UV variability and increasing the efficiency of the simulations, we estimated the relative anomalies *Vi* in *Eery* as a sum of anomalies due to total ozone amount (*v*1*i*,*<sup>j</sup>* (X)) and due to cloud variations (*v*2*i*,*<sup>j</sup>* (*Cl*)). The approach is appropriate, since ozone and cloudiness are located at different levels of the atmosphere. This simple additive method has been used in the UV reconstruction model described in [33], where a good agreement has been shown between results obtained with this method and the observed interannual *Eery* variations. It was also applied for reconstructing the long-term UV variability over Moscow in [10]. The total *Eery* anomalies *Vi* for year *i* can be estimated with the following equation:

$$V\_i = \sum\_j \left( \mathcal{W}\_j(h\_\varepsilon) \left( \nu \mathbf{1}\_{i,j}(X) + \nu \mathbf{2}\_{i,j}(Cl) \right) / \sum\_j \mathcal{W}\_j(h\_\varepsilon) \right) \tag{1}$$

where index *i* and *j* correspond to a year and a month number respectively. *Wj* (*he*) is a weighting function, which is calculated using the effective solar elevation *he*.

For an accurate estimation of the effective solar elevation we used simulations of *h* with 1-hour resolution over a month. An effective solar elevation for each grid was estimated only for hours with *h* > 0. The annual solar elevation changes were accounted in *W* (*he*), using a power law dependence (*Eery* ~ *he* <sup>α</sup>), where α = 2 according to [45]. The same *W* (*he*) correction has been applied to all three datasets (INM-RSHU CCM, ERA-Interim, and TOMS/OMI). In our simulations, we calculated anomalies (*v*1*i*,*<sup>j</sup>* (X) and *v*2*i*,*<sup>j</sup>* (*Cl*)) relative to 2000.

To evaluate the interannual changes in *Eery* due to ozone, we used a well-known Radiation Amplification Factor (*RAF*) technique [46]:

$$\log(Errg\_{i,j}) = -RAF \times \log(\mathbf{X}\_{i,j}) + \mathbf{C}\_{\prime} \tag{2}$$

where *C* is a constant.

In addition, we took into account that the *RAF* for erythemal irradiance depends on solar elevation. We used the following equation according to [47] to parameterize this relationship:

$$RAF(h\_\ell) = -1.10 \times 10^{-4} h\_\varepsilon^{-2} + 1.57 \times 10^{-2} h\_\ell + 0.665\tag{3}$$

As a result, variations *v*1*i*,*<sup>j</sup>* can be written in the following way:

$$w1\_{i,j}(X) = \left(\frac{X\_{i,j} - X\_{2000,j}}{X\_{2000,j}}\right)^{-RAF(h\_\epsilon(j))},\tag{4}$$

where the index 2000 indicates the year 2000, and, hence, *X*2000,*<sup>j</sup>* is a monthly total ozone column in 2000.

In order to characterize both cloud geometry and cloud optical properties effects on *Eery* we applied Cloud Modification Factor for UV spectral region (*CMFUV*), which is often used for this purpose [5,38] and can be estimated from the standard outputs (downward shortwave radiation in clear-sky and in cloudy conditions) of different CCMs and reanalysis datasets. Using these data, one can easily evaluate shortwave radiation cloud modification factor (*CMF*):

$$\mathbb{C}MF = \mathbb{Q}/\mathbb{Q}\_{\mathbb{O}}\tag{5}$$

where *Q* is the surface shortwave radiation in cloudy sky conditions, and *Q*<sup>0</sup> is the surface shortwave radiation in the cloudless atmosphere.

The spectral correction of the cloud modification factor was obtained using accurate model simulations in UV and total shortwave spectrum according to the 8-stream discrete ordinate radiative transfer method in the TUV model [48] and Monte-Carlo simulations [49] for different atmospheric conditions and solar elevations. As a result, we showed that for optically thin cloudiness with *CMF* > 0.95 we can neglect its spectral dependence, while a correction should be used in other conditions. According to [28], this correction is the following:

$$\text{CMF}\_{LV} = (0.1417 \times \sin(h\_{\ell})^2 - 0.175 \times \sin(h\_{\ell}) + 1.054) \text{CMF}^{(-0.04 \times \sin(h\_{\ell})^2 + 0.554 \times \sin(h\_{\ell}) + 0.609)} \tag{6}$$

The proposed method has been successfully tested against the long-term measurements of erythema radiation at the Moscow State University Meteorological Observatory over the 1999–2015 period [28].

To account for total ozone and cloud variations we applied the same method to the different datasets (INM-RSHU CCM, ERA-Interim reanalysis, satellite TOMS/OMI). We should note, that in this study we do not take into account changes in aerosol optical properties and surface reflectivity, which could have some direct and indirect radiative effect on variations in *Eery*, but at much smaller level comparing with ozone and cloudiness over the 1979–2015 period [28].

According to Equations (1) to (6) we evaluated anomalies in *Eery* relative to the year 2000 taking into account either only total ozone or only the cloud modification factor, and considering both of these parameters. After obtaining the *Eery* anomalies, the comparisons between the mean anomaly in *Eery* over the 2000–2015 period and over the 1979–1999 period were performed. For characterizing the statistical significance of the differences in *Eery* between the two periods, we applied a t-test for the means of two independent samples (Welch's *t*-test) at *p* = 95%.

In addition, we performed a trend analysis for *Eery* anomalies using a linear regression model with the least squares approach and estimated its uncertainty using t-test at *p* = 95%.

The details of the UV resources estimations are given in Section 2.4.

#### *2.1. The Description of INM-RSHU Chemical Climate Model*

The INM-RSHU global three-dimensional chemistry–climate model consists of two modules: a dynamical module, which has been developed at the Institute of Numerical Mathematics (INM) of Russian Academy of Science (RAS) [44] and a photochemical module developed at the Russian State Hydrometeorological University (RSHU). Both modules were successfully applied in different international projects on climate and atmospheric gas composition variability [50–52].

The model resolution is 5◦ × 4◦ (longitude x latitude) with 39 vertical levels with variable spacing, up to 0.003 hPa at the model top [43]. The algorithm of the combined model accounts for the interaction between chemical and physical processes at each time step of the model. The chemical module accounts for 74 gas species interacting in 174 chemical reactions, and 46 reactions of photodissociation. A detailed description of the model can be found in [43,44]. For the numerical experiments we used the following datasets: the emissions of the ozone depleting substances were taken from [6]; the change in greenhouse gas emissions was parameterized according to the RCP4.5 scenario from the CMIP5 project; the data on extraterrestrial solar flux variability were taken according to [53], stratospheric aerosol concentrations followed [54]. The simulations were made using different datasets on SST/SIC (Sea Surface Temperature and Sea Ice Coverage, respectively): using the Met Office (as a default variant) [55], the ERA-Interim [42], and the CCM SOCOL (modeling tools for studies of SOlar Climate Ozone Links) datasets [56]. For evaluating *CMFuv*, a spectral correction was applied to the *CMF* data according to Equation (6).

#### *2.2. The Description of ERA-Interim Reanalysis Dataset*

The ERA-Interim reanalysis dataset [42] produced by the European Centre for Medium-Range Weather Forecasts (ECMWF) is a well-known reanalysis dataset available from 1979, continuously updated in real time. The ERA-Interim data-assimilation system utilizes a four-dimensional variational data-assimilation scheme (4D-Var) and a variational bias correction scheme (VarBC) for satellite radiances, that automatically detects and corrects for observation biases [57–60].

The ECMWF ERA-Interim reanalysis dataset, along with data on total ozone content, contains the data on solar radiation for cloudy and clear-sky conditions at the Earth's surface. For calculating *CMF,* we used the following ERA-Interim parameters: daily dose of "surface net solar radiation, clear sky" (NSR, clear sky) and "surface net solar radiation" (NSR), synoptic monthly means. It is easy to show that these parameters can be used for estimating *CMF* using the equation:

$$\text{CMF}\_{i,j} = \frac{\text{NSR}\_{i\bar{j}}}{\text{NSR}\_{\text{clear sky}}\,\,\_{ij}}\,\,\,\,\tag{7}$$

the UV correction was made using the same approach by Equation (6), which has been applied to the *CMF* data obtained from model simulations.

Validation of the ECMWF ERA-Interim ozone dataset has been made for the 1989–2008 period by comparisons with independent ground-based ozone observations [57] and satellite data [58]. In addition to ground-based independent observations, ozone profiles from ozonesondes retrieved from the World Ozone and Ultraviolet Radiation Data Centre (WOUDC) archive were used to validate the ozone vertical distribution. The residuals between the ERA-Interim and the ground-based Dobson total ozone measurements were within ±10% at high latitudes, and within ±5% over other regions with the absence of temporal bias both in the stratosphere and in the troposphere [57]. The validation of ERA-Interim relative to satellite data included total ozone data from TOMS and OMI satellite instruments as well as SBUV and SBUV/2 and ozone profile retrievals from other satellite measurements [58,59]. The ERA-Interim total column ozone uncertainty is typically ±5 DU (about ±2%). As a result, in the latest Scientific Assessment of Ozone Depletion: 2018 [8] on page 3.8 it was stated that "recent reanalysis datasets (Dee et al., 2011; Dragani 2011; and Wargan et al., 2017) have been shown to produce a realistic representation of total ozone".

ERA-Interim global solar irradiances at the ground have been validated with ground-based and satellite data [61–63]. According to a comparison of data from 674 sites, Era-Interim and ground-based measurements were highly correlated (R<sup>2</sup> = 0.97) [61]. Biases relative to observations were smaller compared to other reanalysis data evaluated in this study. According to Table 2 from [61], the mean bias relative to solar radiation measurements from all stations was 11.25 W m−2, and the RMSE (Root-Mean-Square Error) was 27.7 W m−2, while for the best radiation measurements at BSRN (Baseline Surface Radiation Network) sites the mean bias was 4.15 W m−<sup>2</sup> and the RMSE was 19.6 W m<sup>−</sup>2. A similarly small bias of 5.2 W m−<sup>2</sup> was obtained over central Europe [62], while for the French radiation network the mean difference of 9.7 Wm−<sup>2</sup> with *R*<sup>2</sup> = 0.97 was estimated on a daily basis for the 1995–2006 period [63].

Therefore, the estimates of total ozone and solar radiation in the ERA-Interim dataset can be characterized as realistic and suitable for our study.

#### *2.3. Description of the Combined TOMS*/*OMI Satellite Dataset*

The TOMS and OMI satellite datasets were used for comparisons with the *Eery* retrievals from the INM-RSHU CCM and the ERA-Interim reanalysis. For characterizing *Eery* variability, we applied daily erythemal doses (J m<sup>−</sup>2) from the combined satellite dataset. Nimbus7/TOMS and EarthProbe/TOMS Level 3 data were downloaded from https://disc.gsfc.nasa.gov (date of access 01.02.2017) with a spatial resolution of 1◦ × 1.25◦ [64,65]. Level 3 AURA/OMI Erythemal Dose Daily product archives were downloaded from https://giovanni.sci.gsfc.nasa.gov (date of access 16.02.2017) with a spatial resolution of 1◦ × 1◦ [66]. The data were re-gridded to have similar spatial resolution. In addition, we made an updated correction on absorbing aerosol following the methodology described in [17,20,21]. For this purpose, we applied the new aerosol climatology Macv2 for 2005 [41]. The absorbing aerosol optical thickness (τ*abs*) at λ = 320 nm was simulated using the following equation:

$$
\pi\_{\text{abs}} = \tau\_{\text{ext}} (1 - \omega) \tag{8}
$$

where τ*ext* is the extinction aerosol optical thickness; ω is the single scattering albedo.

According to [20,21,67], we evaluated aerosol correction factor *CF*, using the following equation:

$$CF = \frac{1}{1 + K\tau\_{abs}}\tag{9}$$

where the slope *K* weakly depends on solar elevation and aerosol type, and for *h* > 30◦ with typical values of τ*abs* < 0.1 an average value of *K* = 3 can be applied with an error of less than 5% [20].

We should note that in OMI UV retrievals the absorbing aerosol correction had been already applied using the old Macv1 aerosol dataset [20,68]. Hence, for the application of the new Macv2 climatology correction we first removed the previous Macv1 correction from the dataset. For the TOMS

dataset the correction using the new aerosol Macv2 dataset has been applied directly to the standard *Eery* retrievals. We would like to emphasize that both the previous and the new Macv2 corrections do not take into account changes in aerosols over time. The monthly mean τ*abs* distribution with 1◦ × 1◦ grid over the 1980–2015 period have been used as input data for the *CF* estimation. Figure 1 presents relative difference (in %) between the old (Macv1) and the new (Macv2) *CF* aerosol correction factor.

**Figure 1.** Relative difference (in %) between the new (Macv2) and the old (Macv1) aerosol correction factors, *CF* for the central months of the seasons: January (**A**), April (**B**), July (**C**), October (**D**).

Due to smaller absorbing aerosol optical depth over Europe and China in the new aerosol climatology we have relatively higher *CF*'s of up to 10–20% compared with the old Macv1 correction, especially in April and July. The opposite tendency is observed over the Arctic regions due to the increase in τ*abs* in the new climatology, which resulted in 20–30% *CF* decrease almost for all seasons, except January. In other areas there is a slight change in *CF* of about ±5%. Due to changes in absorbing aerosol, there are pronounced differences in UV indices of more than 1 over several regions of Europe, Central Asia and China in July (Figure 2).

**Figure 2.** Absolute difference between the UV indices in 2005 with the application of the updated Macv2 *CF* factor minus UV indices with the standard OMI *CF* factor correction [68] for central months of the seasons: January (**A**), April (**B**), July (**C**), October (**D**). UV index is a widely used unitless quantity, defined by multiplying erythemal irradiance by 40 [1].

In addition, for characterizing the variations in *Eery* due to total ozone changes, the extended observational assimilated database described in [69] has been applied. The assimilated database combines satellite-based ozone measurements from TOMS, GOME, SBUV, and OMI satellite ozone instruments with their corrections from ground-based measurements. Since the development of this database was mainly based on TOMS and OMI total ozone retrievals, we refer to it in this study as the TOMS/OMI ozone dataset. This assimilated database has been also actively used in [8] for ozone trend studies.

We also applied the *CMFuv* climatology obtained previously according to the TOMS Reflectivity measurements (1979–2002) with additional surface albedo correction as described in [14,70]. This dataset has been used for the comparisons with the *CMFuv* climatology evaluated from the INM-RSHU CCM and ERA-Interim datasets. We did not use satellite TOMS/OMI*CMFuv* retrievals for characterizing year-to-year variability, since they have large uncertainty at high latitudes due to high snow surface albedo and low solar elevations [21]. We should mention that the results, obtained using satellite data in other publications [17–19], have been also restricted to 50◦–55◦ latitude to avoid significant uncertainties.

#### *2.4. The UV Resources Simulations*

For characterizing the changes in erythemal UV radiation over the 1979–2015 period and for qualifying the importance of such changes for human health we used the UV resources approach, which was described in details in [14]. As stated above, the UV resources characterize UV radiation from the point of beneficial (vitamin D production) or adverse (sunburn) effects on human health in the categories of UV deficiency, UV optimum and UV excess with the additional account for skin types and open body fraction. The main principle of the approach is in the application of two thresholds. The first threshold is defined by the Minimum Erythemal Dose, *MED*, which depends on skin type

(*k*) according to the Fitzpatrick classification [71]. In this study we consider four skin types, which are typical for Northern Eurasia area. The *MED* for the most UV susceptible skin type 1 is 210 J m<sup>−</sup>2, for skin type 2–250 J m<sup>−</sup>2, for skin 3–350 J m−2, and for skin type 4–450 J m−<sup>2</sup> [71]. For higher UV levels which exceed *MED,* in [14] we proposed to define several subclasses of *UV excess.* They are attributed to the thresholds depending on the WHO (World Health Organization) UV index (*UVI*) categories [72]: *moderate UV excess* category, *high UV excess* category, *very high UV excess* category, and *extremely high UV* excess category, which are respectively related to moderate, high, very high, and extremely high hourly *UVI* (see the details in [14]).

The second threshold is the minimum vitamin D dose production *MvitDD.* A detailed description of the method of its evaluation can be found in [14]. It has been estimated with account for a daily dietary vitamin-D intake of 1000 IU [73] and using an equivalent of vitamin D production of one *MED* of 10,000 IU [74]. Hence, taking into consideration an open body fraction *S*, a relationship between the two thresholds can be written as follows:

$$MvitDD\_k = 0.1 \, MED\_k / S \tag{10}$$

These two thresholds were used for defining different classes of UV resources including the UV deficiency, when no vitamin D is possible to obtain, UV optimum, when it is possible to obtain vitamin D and not possible to get sunburn (erythema), and UV excess, when erythema is obtained. We should note that S was set to 0.25 in this study. This simplified approach provides an opportunity in a first approximation to quantify the UV resources and their spatial and temporal changes.

For evaluating UV resources, absolute *Eery* estimates are needed. For this purpose, we calculated the hourly *Eery* averages centered at noon with 3-min resolution, using the 8-stream discrete ordinate radiative transfer method in the TUV model [48] over the territory of 10.5◦ W–179.5◦ E, 40◦ N–80◦ N with 1◦ × 1◦ increments for year 2000 for clear-sky conditions. In these simulations we took into account total ozone as well as aerosol optical thickness and surface albedo for central days of each month of 2000. To simulate *Eery* in cloudy conditions, we applied the *CMFUV* values. Since *Eery* anomalies in Equation (1) were also simulated relative to 2000, we could easily reconstruct absolute *Eery* changes over the whole period.

Previously, the climatology of UV resources over Northern Eurasia has been obtained in [14]. In this paper we use this approach for evaluating possible changes in spatial distribution of the different UV resources classes since 1979 for different seasons and for various skin types due to ozone and cloud variations.

#### **3. Results**

#### *3.1. Climatologies of Total Ozone and Cloud Modification Factor over Northern Eurasia*

Since total ozone is one of the main factors affecting erythemal UV radiation, it is necessary to calculate both the ozone climatology and its temporal variation with low uncertainty. The comparisons of different ozone climatologies over the 1979–2015 period obtained from INM-RSHU CCM, ERA-Interim Reanalysis and Satellite TOMS/OMI datasets are shown in Figure 3 and Table 1, where the main statistics of total ozone column are presented. The CCM adequately reproduces the absolute values of total ozone amount: the differences between monthly mean CCM total ozone values do not exceed 7% compared with TOMS/OMI ozone retrievals with mean annual difference of about 1% (see Table 1).

**Figure 3.** Climatology of the total ozone content (in Dobson units) over 1979–2015 according to the INM-RSHU CCM, ERA-Interim reanalysis, and TOMS/OMI satellite data. *X*-axis—is longitude, *Y*-axis—is latitude.

**Table 1.** Main statistics of total ozone column X (DU) and cloud modification factor *CMFuv* over Northern Eurasia for the 1979–2015 period \*.


\* For comparison purpose the calculations were made only over the same areas.

There is a good agreement between ERA-Interim and TOMS/OMI mean total ozone datasets with mean differences of 1% or less. There is also a satisfactory agreement in the representation of ozone spatial distribution in different seasons, except October, when a dipole structure of ozone maximum has been obtained with the centers over the Western and Eastern coast of Northern Eurasia, while according to the satellite TOMS/OMI dataset we see only one ozone maximum over Northeastern Asia. There is a high correlation between TOMS/OMI and ERA-Interim Reanalysis datasets since in ERA-Interim ozone is both modeled and assimilated from SBUV, OMI, TOMS, GOME, and SCIAMACHY datasets [57,58].

The effect of cloud on UV radiation is a combination of the effects of cloud optical properties and spatial cloud characteristics (cloud amount). As shown above, cloud modification factor (*CMFuv*) captures both effects simultaneously. Figure 4 shows the spatial distribution of *CMFuv* according to INM-RSHU CCM, ERA\_INTERIM, and TOMS/OMI satellite datasets over the 1979–2015 period. Main statistics of *CMFuv* are shown in Table 1. Noticeable variability in *CMFuv* reflects significant changes of different synoptic processes over Northern Eurasia. In January in midlatitudes cloud modification factors over Europe are much smaller due to the enhanced cyclonic activity (*CMFuv* = 0.3–0.4) and over the central part of the continent the *CMFuv* can reach 0.8–0.9 due to the influence of the Siberian anticyclone, which is characterized mainly by cloudless conditions. In July over the southern European regions one can see high *CMFuv* > 0.9 due to the influence of Azores anticyclone. All datasets are in reasonably agreement, except satellite dataset over Arctic area, where in April and July there is a significant drop in *CMFuv* possibly due to the problems with dividing cloud/albedo effects in the satellite algorithm, that results in unrealistically low *CMFuv* values.

**Figure 4.** Climatology of the cloud modification factor (*CMFuv*) according to different datasets over the 1979–2015 period. Note, that the satellite *CMFuv* climatology has been obtained according to the TOMS data for the 1979–2002 period according to [14]. *X*-axis—is longitude, *Y*-axis—is latitude.

There is a good agreement between all three datasets in July, while in January we observe the *CMFuv* overestimation over the continental areas of Northern Eurasia in the INM–RSHU CCM.

Table 1 also demonstrates the low TOMS/OMI *CMFuv* values compared with the INM–RSHU CCM and ERA-Interim *CMFuv* retrievals, especially in January, in high snow albedo conditions, when the negative bias can reach 20–30%. On average, there is a satisfactory agreement between the mean *CMFuv* values in the INM-RSHU CCM and ERA-Interim datasets with the difference not exceeding 2.5%, except January, when it increased up to 8% (Table 1). We also should mention that previous comparisons with the direct *CMFuv* evaluation from the Moscow State University Meteorological Observatory dataset also demonstrate a satisfactory agreement with the INM–RSHU CCM *CMFuv* values with the exception of winter months [28].

#### *3.2. Changes in Eery Daily Doses due to Ozone and Cloudiness*

#### 3.2.1. Eery due to Total Ozone Variations

There are several factors affecting ozone content in the atmosphere and, hence, noticeably change *Eery* at ground. Due to the Montreal Protocol with its Amendments and Adjustments the emissions of ozone depleting substances (ODS) have been restricted, however, due to significant life time of most of them, their concentrations in the atmosphere are getting lower only in recent years [6,8,10,38]. Along with the anthropogenic effects, the observed variability of some natural factors also could be important in assessing the total variability of ozone. We performed several numerical experiments for evaluating and comparing the role of the main ozone drivers in *Eery* variability over the 1979–2015 period. The detailed description of the numerical experiments is shown in Table 2.



In order to estimate the sensitivity of different factors we simulated the difference in *Eery* between the periods of 2000–2015 and 1979–1999 due to the changes in emissions of ODS, the changes in stratospheric aerosol, in SST/SIC, and solar activity separately and due to their combined effects according to the numerical experiment, when all factors are taken into account simultaneously. The results are shown in Figure 5. Over this period there is a noticeable enhancement in *Eery* due to the ozone loss related to the ODS increase. The largest ozone *Eery* growth of more than 4% is observed over polar regions, which is in accordance to other model simulations [7].

**Figure 5.** Difference in the average of annual Eery anomalies for the period 2000–2015 and the period 1979–1999, due to changes in total ozone column according to effects of different anthropogenic and natural factors. INM-RSHU CCM simulations. Statistically significant difference at *p* = 95% are shown by white hatching. *X*-axis—is longitude, *Y*-axis—is latitude.

Solar activity provides the changes in solar radiation at wavelengths shorter than 242 nm, which influence the ozone photochemical processes. Due to a tendency of decreasing solar activity during the last decades, smaller ozone amount is generated providing small, but positive changes in *Eery* (<1%). The increase in concentration of stratospheric aerosol was mainly observed due to strong volcanic eruptions at the end of the 20th century—El-Chichon (1982) and Pinatubo (1991). They resulted in ozone loss initiated by the heterogeneous processes on the surface of sulphate volcanic aerosol particles in the stratosphere. The absence of intensive volcanic activity in the 21 century provided a small *Eery* decrease during the 2000–2015 period (less than 1%) compared with the 1979–1999.

The changes in SST/SIC have the strongest impact on ozone compared with other natural factors presumably due to their influence on the atmosphere dynamic processes. These changes provide increase in *Eery* up to 2% for annual means in 21th century (2000–2015), especially over the northern Arctic regions. The numerical experiment, which took into account for all factors influencing total ozone column, has revealed a positive change in *Eery* up to 1–2% over northern regions of Eurasia and some small negative effects (less 1%) at the southeastern Pacific Ocean region, which are not statistically significant.

Since the effects of sea surface temperature on ozone are important, we also analyzed the possible changes in *Eery* due to the application of the other SST/SIC datasets in the INM-RSHU CCM (the ERA-Interim [42], and the SOCOL SST/SIC datasets [56]). In these numerical experiments the total ozone column was simulated using the same datasets for ODS, stratospheric aerosol and solar activity as before, but the MetOffice SST/SIC dataset was replaced with ERA-Interim and SOCOL data.

The comparisons of the influence of different SST/SIC on variations in *Eery* have revealed small difference in annual *Eery* behavior, especially between the model runs with the Met Office and the SOCOL SST/SIC. All numerical experiments have a similar tendency of annual *Eery* growth of about 1–2% over the central part of Northern Siberia in 2000–2015 compared with 1979–1999, however, there are large differences in seasonal changes. We would like to emphasize that the influence of SST/SIC on variations in *Eery* should be studied further for understanding the physical mechanism of this phenomenon. We may assume that the changes in SST/SIC may influence the intensity of NAO/AO (North Atlantic Oscillation/Arctic Oscillation) and AAO (Antarctic Oscillation) indices, which is small in zonal averages [75], but according to [8], their variations can explain much of the variability in ozone over local areas [76,77].

In addition, we analyzed variations in *Eery* due to the anthropogenic and natural ozone drivers during the 1979–2015 period using linear trend analysis (Figure 6). Similar to Figure 5, statistically significant positive linear *Eery* trends at *p* = 95% of up to 3%/decade have been obtained according to the numerical experiment with the ODS change (marked as "anthropogenic effect"). The sign of the *Eery* trends due to other natural factors are also similar to *Eery* changes shown in Figure 5.

The numerical experiment, which accounted for all factors, was tested against different datasets. Figure 7 shows the difference in the average of *Eery* anomalies for the period 2000–2015 and the period 1979–1999 due to ozone factor using the *Eery* retrievals according to the INM-RSHU CCM, the ERA-Interim datasets, and TOMS/OMI satellite data. There is a satisfactory agreement between annual results with a 1–2% *Eery* increase due to lower ozone in the 2000–2015 period compared to the 1979–1999 period. The annual Eery changes were not statistically significant, except the *Eery* increase over Eastern Europe and Northeastern Asia according to the ERA-Interim *Eery* retrievals (up to 3%). The discrepancy in *Eery* changes due to ozone in different datasets is much higher for the seasons. While the results for the ERA-Interim and TOMS/OMI satellite datasets are in reasonable agreement for April and July, the agreement for January and October is poor. The *Eery* increase obtained in the INM-RSHU CCM in April is mainly in agreement in sign with ERA-Interim and satellite data with more than 5% increase due to lower ozone content in 2000–2015, however, the spatial pattern of changes in *Eery* does not match other datasets. In July and October large territories are characterized by the decrease in model *Eery* due to more rapid recovery of total ozone compared with that in the ERA-Interim and the TOMS/OMI datasets. Similar tendency in *Eery* change has been also obtained in the detailed analysis of trends in *Eery* over Moscow [28].

**Figure 6.** Decadal trends in annually averaged *Eery* anomalies (%/decade) due to changes in total ozone column according to the effects of anthropogenic (ODS) and different natural factors. INM-RSHU CCM simulations, 1979–2015. Statistically significant trends at 95% are shown by white hatching. *X*-axis—is longitude, *Y*-axis—is latitude.

**Figure 7.** Difference in the average of *Eery* anomalies for the period 2000–2015 and the period 1979–1999 due to total ozone changes calculated from the different datasets. Statistically significant differences at 95% are shown by white hatching. Note that the left-bottom panel is identical to the bottom panel of Figure 5. *X*-axis—is longitude, *Y*-axis—is latitude.

In addition, Figure 8 presents time series of *Eery* anomalies at 48◦ N for longitudes between 0◦ E and 150◦ E based on the INM-RSHU CCM, ERA-Interim, and TOMS/OMI satellite datasets. The results from the three datasets agree within a few percent and suggest that *Eery* has increased between 1979 and 1995. In contrast, *Eery* has not changed perceivably between 2000 and 2015. These results corroborate the conclusion from Figure 8 that *Eery* was lower during 1979–1999 compared to 2000–2015.

**Figure 8.** The time series of the *Eery* relative anomalies against 2000 at 48◦ N from 0◦ E to 150◦ E according to the INM-RSHU CCM, ERA-Interim, and TOMS/OMI satellite datasets.

Figure 9 shows the spatial distribution of the decadal trend in *Eery* caused by changes in total ozone over the period 1979–2015, calculated from the INM-RSHU CCM, ERA-Interim, and TOMS/OMI satellite datasets. For annual means (last row in Figure 9), trends computed from the three datasets are generally consistent and range between −0.5% and 2% per decade. However, similar to the results shown in Figure 7, the discrepancy is much larger for the seasonal datasets. In April a pronounced increase in *Eery* (up to 3% per decade) is observed according to all datasets, which is statistically significant over vast areas. Results for the ERA-Interim and TOMS/OMI satellite datasets agree again reasonably well, while results from the INM-RSHU CCM model in most cases do not capture the spatial patterns of the observations.

#### 3.2.2. *Eery* due to Cloud Variations

We analyzed *CMFuv* variations retrieved from the INM-RSHU CCM numerical experiments as well as from the ERA-Interim dataset. The seasonal and annual differences in the *CMFuv* values between the 2000–2015 and 1979–1999 periods are shown in Figure 10. In general, the *CMFuv* differences vary between these periods within ±5%. However, in April over Central and Eastern Europe and in July over the northern Atlantic, the European territory of Russia, Central Siberia and Northeastern Asia the difference is higher than 5% and sometimes exceeds 10% according to the ERA-Interim dataset. These estimates are in a good agreement with the *CMFuv* trends over Moscow [78]. Model simulations provide much smaller changes in *CMFuv*. Similar tendencies are observed only over the Atlantic and Northeastern Asia, but at much smaller level. In contrast, in July over the Northern Arctic region there is a statistically significant *CMFuv* decrease, which has been obtained both according to the INM-RSHU CCM and the ERA-Interim datasets. This may happen due to increasing temperature in this northern area as a result of global warming and intensification of the cyclonic processes and, hence, increase in water vapor content. We should emphasize that these results are also in agreement with the independent model simulations [11]. In addition, in [79] it was mentioned that by the end of the 21 century "cloud cover will increase at high latitudes by up to 5% but will decrease at low latitudes (<~30◦) by up to 3%". According to [80] a UV reduction of 10–15% is projected by 2100 due to increases

in cloudiness over some northern high-latitude regions and over Antarctica. Hence, the tendency in increasing of cloud cover in summer Arctic conditions is successfully retrieved from INM-RSHU CCM. At the same time, the INM-RSHU CCM does not reproduce well positive changes in *CMFuv* over midlatitudes. For annual data, one can see even a prevailing small model *CMFuv* decrease compared with the *CMFuv* increase of up to 5% over vast midlatitude areas according to the ERA-Interim dataset.

**Figure 9.** Spatial distribution of the decadal trend in *Eery* caused by changes in total ozone over the period 1979–2015, calculated from the INM-RSHU CCM, ERA-Interim, and TOMS/OMI satellite datasets. Statistically significant differences at 95% are shown by white hatching. Note that the left-bottom panel is identical to bottom panel of Figure 6. *X*-axis—is longitude, *Y*-axis—is latitude.

The additional linear trend analysis also demonstrates high *CMFuv* increase (up to +4–6% per decade) over several areas (East-European Plain, Northeastern Asia, and local areas in Siberia) according to the ERA-Interim dataset (Figure 11), while model positive linear decadal trends are much smaller. Over the Arctic basin negative decadal trends of about 2–4%/decade have been revealed from both model and ERA-Interim *CMFuv* retrievals in April and July. This is also in agreement with the results obtained in [11]. We should also note that the model does not reproduce positive trends in *CMFuv* over Central Europe, Central Siberia and Northeastern Asia in July, and over the Northeastern Arctic regions in October.

**Figure 10.** Difference in the average of *CMFuv* anomalies for the period 2000–2015 and the period 1979–1999 according to the INM-RSHU CCM simulations and the retrievals from the ERA-Interim dataset. Statistically significant differences at *p* = 95% are shown by white hatching. *X*-axis—is longitude, *Y*-axis—is latitude.

**Figure 11.** The decadal trends in *Eery* due to changes in cloudiness according to the INM-RSHU CCM simulations and the retrievals from the ERA-Interim dataset over the 1979–2015 period. Statistically significant trends at 95% level are shown by hatching. *X*-axis—is longitude, *Y*-axis—is latitude.

#### 3.2.3. *Eery* Changes due to Joint Influence of Total Ozone and Cloud Variations

According to the method, described in Section 2 (Equations (1) to (6)), we simulated variations in *Eery* due to both factors: ozone and cloudiness. We applied this approach for both ERA-Interim and INM-RSHU CCM datasets. In addition, we made the comparisons with TOMS/OMI daily erythemal doses dataset with the new Macv2 aerosol correction. Figure 12 shows seasonal and annual differences in *Eery* anomalies between the 2000–2015 and 1979–1999 periods due to the combined effects of total ozone and cloudiness. For annual means (last row in Figure 12), results for the CCM model show increases by 0–3% in 2000–2015, while changes for the ERA-Interim dataset range between −2% and 5%, with few exceptions. In July the spatial changes in *Eery* from the Era-Interim dataset exhibit a quasi-wave spatial structure in midlatitudes with a significant increase of about 10–15% over Europe, Central Siberia, and Northeastern Asia. A similar structure, albeit less pronounced, is also apparent in satellite retrievals, and the CCM model. The agreement in the spatial patterns of the datasets is generally poor for other months, except areas over oceans. This disagreement may partly be caused by the interaction between clouds and high surface albedo from snow cover (e.g., see discussion in [67]). This interaction is particularly a problem over northern areas that are affected by snow cover from October through June.

**Figure 12.** Seasonal and annual differences in *Eery* anomalies between the 2000–2015 and 1979–1999 periods due to combined effects of changes in total ozone column and cloudiness according to the INM-RSHU CCM model, the ERA-Interim dataset, and Eery retrievals from the TOMS/OMI satellite dataset with the additional Macv2 aerosol correction. Statistically significant difference at 95% are shown by hatching. Note, that we do not present UV retrievals from the TOMS/OMI dataset for January and for annual mean due to the problems in UV retrievals in conditions with high surface snow/ice albedo. *X*-axis—is longitude, *Y*-axis—is latitude.

Figure 13 presents *Eery* decadal trends due to the combined effect of total ozone and cloudiness according to the INM-RSHU CCM, ERA-Interim simulations, and *Eery* retrievals from the TOMS/OMI satellite dataset with the additional aerosol Macv2 correction. For annual means (last row in Figure 13), trends calculated from the CCM and ERA-Interim datasets generally agree, range between 0% and 3% per decade, and are statistically significant over large areas. The noticeable positive trends in *Eery* of up to 8% per decade in July, and up to 5% per decade in April are observed in the ERA-Interim dataset over several areas in Eastern Europe, Central Siberia, and Northeastern Asia. In July, similar trends in *Eery* of up to 8% per decade can be found according to the satellite data over the Pacific Ocean. We also see a better agreement in July, when the spatial quasi-wave structure of positive *Eery* trends in midlatitudes from the ERA-Interim is similar to that obtained from the satellite data. Spatial patterns

of the CCM and ERA-Interim results agree qualitatively, but trends calculated with the CCM model are generally smaller.

**Figure 13.** Spatial distribution of the decadal trend in *Eery* caused by the combined effect of total ozone and cloudiness calculated from the INM-RSHU CCM, ERA-Interim, and TOMS/OMI satellite datasets. Statistically significant differences at 95% are shown by hatching. Note, that we do not present UV retrievals from the TOMS/OMI dataset for January and for annual mean due to the problems in UV retrievals in conditions with high surface snow/ice albedo. *X*-axis—is longitude, *Y*-axis—is latitude.

#### *3.3. Changes in UV Resources due to the Changes in Ozone and Cloudiness over Northern Eurasia*

Using the method described in Section 2.4, we simulated the UV resources for 1979 and 2015 and their changes over this period according to the retrievals from the ERA-Interim data for different skin types. The *Eery* retrievals from the ERA-Interim dataset were chosen for the following reasons. For estimating *Eery*, we need the reliable products for both total ozone and cloud modification factor. A large territory of Northern Eurasia is located in high latitudes with small solar elevations and snow cover dominating not only in winter, but during spring and fall, when satellite *Eery* retrievals cannot be used. In Section 2.2 the Era-Interim total ozone and downward shortwave radiation products were shown to have a good quality [57–63]. Therefore, we simulated UV resources using the *Eery* retrievals from this dataset. To avoid local features in *Eery* variations in 1979 and 2015, we applied a linear regression model to characterize the *Eery* changes from 1979 to 2015.

Figure 14 presents the UV resources distribution over Northern Eurasia for 1979 and 2015 and their difference for skin types 1–4. According to [72], skin type 1 is most susceptible to erythema, while skin type 4 is characterized by small sun-sensitivity. We revealed that for the all skin types due to the reduction of ozone and clouds, there is a noticeable geographical shift of UV categories over large areas, especially in April and July. A noticeable reduction of the UV optimum area and its replacement by the moderate UV excess conditions is observed in April for skin types 1 and 2. There is also a shift from moderate to high UV excess conditions at lower latitudes for these skin types. However, for skin types 3 and 4 in April we see favorable changes from UV deficiency to the UV optimum conditions over polar regions, and adverse change of UV optimum to UV excess conditions at the south (see Figure 14). In July in central Arctic region there is, on the contrary, a change from UV excess to UV optimum conditions for skin types 2 and 3, due to an increase in cloudiness. In October favorable changes are observed for skin types 3 and 4 with UV optimum area shift towards the north, while for skin types 1 and 2 there is the replacement of UV optimum by moderate UV excess conditions at the south.

Skin type 1

#### Skin type 3

**Figure 14.** UV resources for 1979 and 2015 and changes in their distribution between 2005 and 1979 according to the retrievals from the ERA-Interim data for different skin types. *X*-axis—is longitude, *Y*-axis—is latitude.

#### **4. Discussion**

In this paper we evaluated the separate effects of total ozone and cloudiness and their combined effect on temporal variability of UV erythemal daily doses and UV resources according to the *Eery* retrievals from the ERA-Interim, the INM-RSHU CCM datasets, and TOMS/OMI satellite measurements with the Macv2 aerosol correction over Northern Eurasia for the 1979–2015 period.

The numerical experiments using INM-RSHU CCM model allowed us to estimate the role of main ozone drivers (ODS, stratospheric aerosol, sea surface temperature and ice coverage, solar activity) and to evaluate their impact on *Eery*. We showed that the most important factor affecting *Eery* is the concentration of ODS, the increase of which led to the significant ozone loss. Among natural factors, stratospheric aerosol and SST/SIC are the most important ones. The application of the Met Office and the SOCOL SST datasets in the INM-RSHU CCM provided the similar annual *Eery* increase in 2000–2015 over Central part of Northern Siberia, however, further studies should be conducted for understanding the physical mechanism of this phenomenon. On the whole, there is a positive change in the modeled *Eery* of up to 1–2% at the northern regions of Eurasia in 2000–2015 compared with the 1979–1999 period. Maximum positive linear *Eery* trends (up to 3% per decade) due to ozone loss were observed in April.

Linear *Eery* trends due to changes in cloudiness were much more pronounced and reached 4–8% per decade over some areas (East-European Plain, Northeastern Asia and local spots in Siberia) according to the ERA-Interim dataset. Both INM-RSHU CCM and ERA-Interim *Eery* retrievals have revealed negative *Eery* trends due to cloudiness of about 2–4% per decade in the central Arctic basin in April and July. This decrease in *CMFuv* during the last decades is connected with warming of the lower troposphere and large reduction of sea-ice-sheet, increasing evaporation and water vapor content [81]. Similar changes were obtained in [11], where it was mentioned that "the primary drivers of these changes are increasing concentrations of greenhouse gases (GHGs) and, for the southern hemisphere, the Antarctic ozone 'hole'". It should be emphasized that according to the model estimations *CMFuv* changes were much smaller than those obtained from ERA-Interim and satellite data. According to the ERA-Interim *Eery* retrievals, the joint ozone and cloud effect provided up to 6–8% increase in *Eery* per decade in July and April over several areas in Eastern Europe, Central Siberia, and Northeastern Asia. Annual *Eery* linear trends comprised about 2% per decade according to the ERA-Interim dataset and were statistically significant over large areas.

The analysis of the UV resources for 1979 and 2015 and their changes over this period showed a noticeable geographical shift of UV categories, especially in spring and summer over large areas due to the reduction of ozone and clouds. For skin types 1 and 2, which are most susceptible to sunburn, favorable UV optimum conditions were substituted by UV excess conditions over large areas. At the same time in April favorable UV optimum instead of UV deficiency conditions were observed for skin types 3 and 4 in the north but in the south the areas covered by the UV optimum have been replaced by the UV excess conditions. In some Arctic regions, in July, there was a change from UV excess to UV optimum for skin types 2 and 3, due to the increase in cloudiness.

In our assessment of UV resources, we did not take into account changes in the surface area of the skin that is exposed to sunlight or possible changes in the behavior of people, which could result from climate change and increasing temperatures. Instead, we focused on variations in UV resources due to the changes of the most important geophysical factors: total ozone and cloudiness.

#### **5. Conclusions**

We analyzed temporal variability of UV erythemal daily doses and UV resources due to ozone and cloudiness, according to different datasets, using *Eery* retrievals from the ERA-Interim reanalysis, INM-RSHU CCM, and the TOMS/OMI satellite measurements, with the updated aerosol correction over Northern Eurasia for the 1979–2015 period.

We showed that, according to all datasets for spring and summer clear-sky conditions, there was a pronounced trend in *Eery* of up to +3% per decade due to ozone loss, which was statistically significant over large areas in Northern Eurasia.

The INM-RSHU model experiments have confirmed that the largest effect on ozone and *Eery* stem from the largest impact of anthropogenic emissions of ozone-depleting substances. Additional factors include volcanic aerosol and the SST/SIC on ozone and, hence, on *Eery* changes. The utilization of the different SST/SIC datasets in the INM-RSHU CCM showed similar annual mean *Eery* increase over the polar region in Siberia.

The INM-RSHU CCM, as many other CCM models (see, for example, discussion in [37]), did not reproduce an observed significant positive change in *CMFuv* during the last decades, which can reach up to 4–8% per decade, according to the ERA-Interim dataset over several areas in Northern Eurasia.

According to the *Eery* retrievals from the ERA-Interim dataset positive *Eery* trends reached 6–8% per decade over Eastern Europe, several regions in Siberia and Northeastern Asia due to the joint impact of ozone and cloudiness. In the central Arctic region, negative *Eery* trends were observed in summer due to *CMFuv* decrease, which is in agreement with findings in other studies [11].

The simulations of changes in UV resources from 1979 to 2015 have shown that there is a shift toward higher UV categories, especially in April and July, for skin types 1 and 2, which are most susceptible to sunburn. In April, this phenomenon was expressed in a noticeable reduction of the UV optimum area and its replacement by the moderate UV excess conditions. In addition, a shift from moderate to high UV excess conditions is observed. However, during summer months in the central Arctic region, there is a change from UV excess to UV optimum conditions for the skin types II and III due to a decrease in *CMFuv*.

This work is a part of an on-going project, within which we plan to perform additional numerical experiments with the INM-RSHU-CCM taking into the account temporal changes in aerosol and surface reflectivity as well as aerosol–cloud interaction, which possibly help in evaluating real *CMFuv* long-term changes over 1979–2015 period.

**Author Contributions:** The paper was initiated and written by N.E.C. All other authors had valuable contribution in visualization of the results and data analysis (A.S.P., E.Y.Z., and E.V.V.) and in providing the data of model experiments with the INM-RSHU CCM (S.P.S. and V.Y.G.); E.V.V. and E.Y.Z. helped with the technical preparation of the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was partly supported in part by the Russian Foundation for Basic Research, grant no. #18-05-00700, and state budget program AAAA-A16-116032810086-4.

**Acknowledgments:** We would like to thank A. Arola for providing the aerosol dataset, which was used in the standard OMI *Eery* retrievals. We would also like to thank Bodeker Scientific, funded by the New Zealand Deep South National Science Challenge, for providing the combined NIWA-BS total column ozone database. In addition, the useful comments of the anonymous reviewers are highly acknowledged.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Discontinuities in the Ozone Concentration Time Series from MERRA 2 Reanalysis**

#### **Peter Krizan \*, Michal Kozubek and Jan Lastovicka**

Institute of Atmospheric Physics, Czech Academy of Sciences, 14100 Prague, Czech Republic; kom@ufa.cas.cz (M.K.); jla@ufa.cas.cz (J.L.)

**\*** Correspondence: krizan@ufa.cas.cz

Received: 17 October 2019; Accepted: 12 December 2019; Published: 14 December 2019

**Abstract:** Artificial discontinuities in time series are a great problem for trend analysis because they influence the values of the trend and its significance. The aim of this paper is to investigate their occurrence in the Modern-Era Retrospective analysis for Research and Applications, version 2 (MERRA 2) ozone concentration data. It is the first step toward the utilization of the MERRA 2 ozone data for trend analysis. We use the Pettitt homogeneity test to search for discontinuities in the ozone time series. We showed the data above 4 hPa are not suitable for trend analyses due to the unrealistic patterns in an average ozone concentration and due to the frequent occurrence of significant discontinuities. Below this layer in the stratosphere, their number is much smaller, and mostly, they are insignificant, and the patterns of the average ozone concentration are explainable. In the troposphere, the number of discontinuities increases, but they are insignificant. The transition from Solar Backscatter Ultraviolet Radiometer (SBUV) to Earth Observing System (EOS) Aura data in 2004 is visible only above 1 hPa, where the data are not suitable for trend analyses due to other reasons. We can conclude the MERRA 2 ozone concentration data can be used in trend analysis with caution only below 4 hPa.

**Keywords:** merra ozone data; discontinuities in reanalysis time series; trend analyses

#### **1. Introduction**

Ozone is a very important trace gas in the atmosphere because it protects the earth's biota from harmful ultraviolet radiation. The high scientific interest in atmospheric ozone arose near the beginning of the 1970s in connection with the problem of the possible supersonic transport impact on the ozone layer. After discovering the Antarctic ozone hole [1], the ozone became an important topic of atmospheric research. The origin of the ozone hole is in the chemical reactions of anthropogenic halogen radicals, which destroy ozone [2]. These results led to signing the Montreal Protocol and its amendments, which stopped the production of the ozone-depleting substances (ODS). This protocol has a positive influence on the ozone layer due to the decrease in ODS [3]. There are some signs of the ozone recovery, especially in the upper stratosphere [4]. In addition, all model studies predict future recovery of the ozone layer [5]. The concentration of ODS is not the only factor that influences the ozone layer. In addition, the greenhouse cooling of the stratosphere [6] and changes in the Brewer–Dobson circulation [7] have an impact on the ozone concentration. In such situations, proper trend analysis is necessary for the understanding of the ozone layer behaviour. Trend analysis was done from ground-based data [8,9]. These data are measured at one point, and the number of ground-based measurements of ozone is insufficient, especially in the Southern Hemisphere. Satellite ozone data have broader coverage, and these data are widely used in trend analysis in ozone research [10]. But in some areas, it is impossible to measure ozone (polar night, below dense clouds) by satellite.

On the other hand, the ozone data from reanalysis are temporally and spatially homogeneous, but there is a big question of the suitability of these data for trend analysis due to the occurrence

of discontinuities in them [11]. They can be caused by satellite or instrument replacement or by assimilation of not homogenous basic parameters. [12] considered the evaluation of trends from reanalysis to be a challenge. The aim of this paper is to detect grid points in which discontinuities occur with the help of the Pettitt homogeneity test in the MERRA 2 ozone data in the period 1980–2017. Our study can help us to identify where it is possible to use the MERRA 2 reanalysis for trend study. This paper is divided into the following sections: Section 2 describes the data and method, Section 3 provides results, Section 4 discusses the results, and Section 5 concludes our results.

#### **2. Data and Method**

In this paper, we used MERRA 2 ozone monthly mean data from 500 hPa up to the top of the model (0.01 hPa) in the period 1980–2017. MERRA 2 reanalysis covers the whole satellite era. It includes substantial upgrades and changes to the data assimilation system and input data in comparison to previous MERRA reanalysis. New constraints are applied to ensure the conservation of global dry-air mass and to close the balance between surface water fluxes (precipitation minus evaporation) and changes in total atmospheric water [13]. The modified gravity wave scheme substantially improves the model representation of the quasi biennial oscillation (QBO) [14,15]. The assimilation of Microwave Limb Sounder (MLS) temperature retrievals at high-pressure levels (higher than 5 hPa) should improve the reanalysis at upper levels. The assimilation of MLS stratospheric ozone profiles and Ozone Monitoring Instrument (OMI) column ozone since the beginning of the Aura mission in late 2004 also improve the representation of fine-scale ozone features, especially in the region around the tropopause [16]. MERRA and MERRA 2 use the Three Dimensional Variational (3D VAR) assimilation process. MERRA 2 uses regular latitude–longitude grids from 1000 to 0.01 hPa (1/2◦ latitude × 5/8◦ longitude). Detail description of all reanalyses (included MERRA 2) can be found in ref [17,18].

We did our analysis for each grid point, we did not use any spatial (longitudinal and latitudinal) averages, because, during every averaging, some information is lost. In each grid and each layer, we used the time series of the ozone concentration and applied the Pettitt homogeneity test [19] to look for discontinuity in it. In each grid, the Pettitt test estimates only one main discontinuity, so this procedure is not able to detect multiple discontinuities. The Pettitt test is a nonparametric test. It has been developed to test homogeneity against the shift in a time series. In this approach, a shift point is at T; indicates that the time series can be divided into two subsequences xt (t = 1..T) and xt (t = T + 1,.. N). Thus, probability distribution functions F1(X) and F2(X) can be associated with the two subsequences. In practice, the null hypothesis H0 is F1(X) is equal to F2(X), and the alternative hypothesis H1 is F1(X) is not equal to F2(X). This method also involves a comparison of the observations so that:

Di,j = sgn(xi − xj), Di, j = 1 if xi > xj, Di, j = 0 if xi = xj, Di, j = −1 if xi < xj. We must compute the statistic: UT, N = t = 1..T <sup>t</sup> <sup>=</sup> T, N D*i*,*j*, t = 1..N

Using the theory on statistic ranks, another KN variable is derived from UT, N This new variable is defined:

KN max abs(UT,N) T = 1..N − 1.

For the test, a probability of exceeding is fixed for a threshold value k given by the formula: P (KN>K) -6exp(−6k2)/(N3 <sup>+</sup> <sup>N</sup>2).

The null hypothesis, H0 is rejected if the probability of exceedance given in P is less than the significance level α for a one-sided statistic test.

This test is widely used in the climatological research especially for precipitation and temperature analysis [20–22]. We searched for temporal occurrence of discontinuities, so in each grid with a discontinuity, we looked for the year of its occurrence, and we were able to detect the discontinuity occurrence in a certain year. We are interested in the spatial distribution of discontinuities as well. So in each layer and month, we constructed the map of their occurrence.

The Pettitt test tells us the year in which the discontinuity in time series occurs. But it does not say how big the discontinuity is or how it can affect trends. Small discontinuities have little impact on the trend analyses. On the other hand, the large one can have a strong trend impact, so we must divide the discontinuities according to their size. We tried to identify which ones were significant (in our case big enough to impact the trend) or insignificant according to this rule: Suppose we have a time series with the length L, and in year x, we observe the discontinuity. We compute the difference between the average before the year x and after this year. If this difference is larger than the variance of time series, we can say that the discontinuity is significant and will have some impact on trends, and we should be careful using this grid point in a trend analysis. This rule should be used mainly in areas with big variance because if the variance is very small, even the small discontinuities can be identified as significant. Another possibility for assessing the size of the discontinuities is the use of the Mann–Whitney test [23] for detecting significant ones. It would be very interesting to compare the results of both methods, and this topic could be solved in the future paper.

#### **3. Results**

#### *3.1. Patterns of Ozone Concentration*

In Figure 1, one can see the map of ozone concentration at 0.1 hPa (upper panel) and at 0.5 hPa (lower panel). In the uppermost model layer, we see the strange pattern of the ozone concentration, which cannot be considered as real. At 0.5 hPa, this strange pattern was also present, but we also see the maximum in ozone concentration above the Aleutian Islands. At 1 hPa in January (Figure 2 upper panel), the strange unrealistic pattern ceased, but the maximum in the ozone concentration appeared over the Aleutian Islands. In July (Figure 2 lower panel), this maximum was situated in the subpolar latitudes of the Southern Hemisphere. At 2 hPa, this maximum was weaker and at 3 hPa disappeared. In the upper stratosphere (Figure 3 upper panel), we observed the equatorial maximum and the polar minimum of ozone concentration. In this layer, the solar radiation is the main factor that drives the ozone concentration. In the lower stratosphere (Figure 3 lower panel), the maximal concentration of ozone was seen in the high latitudes and the minimal one over the equator, which is given by the Brewer–Dobson circulation. The lower stratospheric maximum was persistent in the upper troposphere (Figure 4 upper panel). At 500 hPa (Figure 4 lower panel), this tropospheric maximum disappeared, and we observed maximal ozone concentration over the polluted areas of China and the United States of America. So, we found several ozone concentration patterns from 500 hPa up to the top of the model: tropospheric pattern, lower stratospheric, middle stratospheric pattern, maximum pattern, and strange pattern in the uppermost model layers. The pressure ranges in which these patterns occurred are given in Table 1. The strange pattern occurred above 0.5 hPa. The pattern with a maximum above the subpolar latitudes was seen from 0.7 hPa down to 2 (3) hPa. From October to March, this maximum was situated in the Northern Hemisphere, while from April to October in the Southern Hemisphere. In October, both maxima were present, but the northern maximum was much stronger than the southern one. Below this pattern, the upper stratospheric one with the maximum over the equator and the minimum over high latitudes of both hemispheres could be seen. The lower limit of this pattern was about 20 hPa in each month, and the upper limit was monthly dependant (from 7 to 3 hPa). The upper limit for the lower stratospheric pattern was 30 (40) hPa in each month, while this pattern had its lower edge deep in the troposphere (Figure 4 upper panel) in the Northern Hemisphere and from December to June also in the Southern Hemisphere. From July to November, we saw a secondary ozone concentration maximum also in the subpolar latitudes of the Southern Hemisphere, which is given by the fact the ozone hole starts to form in polar latitudes and so this maximum appeared as a consequence of this formation. If the ozone-destroying reactions

did not act, one could expect a higher concentration of ozone in the polar latitudes of the Southern Hemisphere, and so no subpolar maximum occurs.

**Figure 1.** Average ozone concentration (kg/kg) in January at 0.1 hPa (**upper panel**) and at 0.5 hPa (**lower panel**).

**Figure 2.** The same as Figure 1, but for January at 1 hPa (**upper panel**) and for July at 1 hPa (**lower panel**).

**Figure 3.** Average ozone concentration (kg/kg) in January at 10 hPa (**upper panel**) and at 50 hPa (**lower panel**).

**Figure 4.** Average ozone concentration (kg/kg) in January at 300 hPa (**upper panel**) and at 500 hPa (**lower panel**).

**Table 1.** The vertical extent of ozone concentration pattern in a certain month (satellite—satellite-like pattern), max—the pattern with maximum over the Aleutian Islands (Al) or over Antarctica (An), UST—upper stratospheric pattern (maximum over equator and minimum over polar latitudes), LST—lower stratospheric pattern (maximum over the polar latitudes and minimum over equator), SH max—another maximum over the subpolar latitudes of the Southern Hemisphere.


#### *3.2. Temporal Occurrence of Discontinuities*

Now we deal with the temporal occurrence of discontinuities in the MERRA 2 ozone data. When we look at Table 2, we see in the uppermost model layer in each month the unimodal temporal distribution with a maximum in 1993 (U93 Figure 5—upper panel for January). Going down this unimodal distribution change, the bimodal one with the main maximum was observed in 1993, and the secondary maximum was observed in 2004(B**<sup>93</sup>**,04). In lower layers, the bimodal distribution with the main maximum in 2004 and the secondary one in 1993 (B93,**<sup>04</sup>**) could be seen. This transformation ended at 1 hPa (Figure 5 lower panel for January), where we observed unimodal distribution again, but the maximum occurred in 2004 in each month. Below 2 hPa down to 200 hPa, we observed flat temporal distribution of discontinuities (F) with no sharp maximum (Figure 6 upper panel for January) in each month. In January, the upper distribution border lay in 5 hPa. Below 200 hPa, we again observed unimodal temporal distribution with a maximum about the early 1990s (U93T, Figure 6—lower panel for January), but this distribution was broader than distribution U<sup>93</sup> in the upper

model layers. When we observed the unimodal distribution of all discontinuities, the distribution of the significant ones must be similar. In the case of bimodal distribution, the situation was different. We could observe the bimodal distribution of all discontinuities, but the distribution of their significant counterparts can be unimodal. This was the case of bimodal distributions from 0.4 hPa and 0.5 hPa when we observed the bimodal distribution of all discontinuities, but the 1993 maximum was formed predominantly by the insignificant ones, and the 2004 maximum consisted predominantly of the significant discontinuities (Figure 7). Thus the distribution of the significant discontinuities in these cases was bimodal, but the maximum in 2004 was much stronger than that in 1993. Similarly, in the troposphere, the distribution U93T was seen only for all discontinuities, but for the significant ones, we observed nearly flat distribution due to the fact the 1993 maximum consisted predominantly of the insignificant discontinuities (Figure 8).

**Table 2.** Temporal distribution of discontinuities at certain layers and certain months (U93—unimodal distribution with maximum in 1993, B93,04—bimodal distribution with two maxima in 1993 and 2004, the height of maxima was comparable, B**<sup>93</sup>**,04—bimodal distribution, maximum in 1993 was higher than that in 2004, B93**,04**—bimodal distribution, maximum in 1993 was lower than that in 2004, U04—unimodal distribution with maximum in 2004, F—flat distribution with no maxima, U93T—unimodal distribution with wide maximum in 1993 which occurred in the troposphere).


**Figure 5.** Temporal distribution of all discontinuities in January at 0.1 hPa (**upper panel**) and at 1 hPa (**lower panel**).

**Figure 6.** The same as Figure 5, but for 50 hPa (**upper panel**) and for 300 hPa (**lower panel**).

**Figure 7.** Temporal distribution at 0.5 hPa of all (**upper panel**) and significant (**lower panel**) discontinuities.

**Figure 8.** The same as Figure 7, but for 300 hPa.

#### *3.3. Vertical Profile of the Discontinuities Occurrence*

In Figure 9 (upper panel), we see the vertical profile of the percentage occurrence of all discontinuities. One hundred percent means that the discontinuity occurred in each grid point (207,936 points), and 0% means there were no discontinuities in time series over the globe. The vertical profile pattern was similar in each month. In the uppermost model layer at each grid, the discontinuity occurred. Going down, the share of all discontinuities (Figure 9 upper panel) declined to 0.7 hPa. In 1 hPa, we observed their local maximum. Below this layer, the number of discontinuities declined to about 10 hPa, where we observed local occurrence minimum. From 10 hPa down to 200 hPa, the discontinuities occurred in 40%–60% of all grids. In the troposphere, the share of all discontinuities increased.

**Figure 9.** Vertical profiles of percentage occurrence of all (**upper panel**) and significant discontinuities (**lower panel**) in all months.

The vertical profiles of the significant discontinuities were similar (Figure 9 lower panel), but there were some interesting features in them. The local minimum about 0.7 hPa was much deeper than in the case of all discontinuities. The number of significant discontinuities was nearly zero about 10 hPa, and the tropospheric maximum, which was seen in the case of all discontinuities, was not present for the significant ones. With the increasing difference between the percentage of all discontinuities and the percentage of the significant ones, the share of significant discontinuities decreased.

#### *3.4. Spatial Occurrence of Discontinuities*

Now we deal with the spatial occurrence of discontinuities over the globe. We present the results for January, but the results for other months were similar. Figure 10 (upper panel) shows the geographical distribution of all discontinuities at 0.1 hPa, and in the lower panel, the distribution of the significant ones is displayed. In this layer, there were discontinuities in all grids (orange colour in the upper panel of Figure 10), and their reduction for the significant ones was small. At 3hPa (Figure 11), the situation was different. The number of all discontinuities (upper panel) was smaller, and the number of their significant counterparts (lower panel) was strongly reduced. From 5 hPa down to 250 hPa, the number of significant and insignificant discontinuities were relatively low, and at these layers, almost there were no significant discontinuities (Figures 12 and 13). From 300 hPa down to 500 hPa (Figure 14), the number of all discontinuities was very high, but the number of significant ones was low.

**Figure 10.** Geographical distribution of all (**upper panel**) and the significant (**lower panel**) discontinuities for January at 0.1 hPa (red—discontinuities, yellow—no discontinuities; upper panel all points have discontinuities.).

**Figure 11.** The same as Figure 10, but for 3 hPa.

**Figure 12.** Geographical distribution of all (**upper panel**) and the significant (**lower panel**) discontinuities for January at 10 hPa.

**Figure 13.** The same as Figure 12, but for 250 hPa.

**Figure 14.** Geographical distribution of all (**upper panel**) and the significant (**lower panel**) discontinuities for January at 500 hPa.

Now we search geographical areas where the discontinuities occur more often. This analysis will be done separately for all and significant discontinuities. At each layer and each grid, we computed for all months the sum of cases when we observed discontinuity in a certain grid. When the sum equals 12, it means that in this grid, the discontinuity was observed in each month. On the other hand, when the sum equals 0, it means there was no discontinuity in this grid in every month. The higher the sum, the more frequent occurrence of discontinuities was seen in this grid. At the uppermost model layers, the spatial distribution of the sum was layer-dependant (from 0.1 hPa down to 1 hPa). This can be illustrated in Figures 15 and 16, where in the case of the insignificant discontinuities (lower panel) we observed at 0.3 hPa, a higher sum in the tropics than in the middle latitudes, but at 0.4 hPa (one layer below), the opposite was true. We observed a decrease in the sum from 1 hPa (Figure 17) down to 5 hPa (Figure 18) due to lower numbers of discontinuities. From 5 hPa down to 250 hPa, this sum for all and the significant discontinuities was relatively low, because the occurrence of discontinuities was rare. (Figure 19). In the troposphere (300 hPa Figure 20), we observed the increase in this sum, especially over the Pacific, in the case of all discontinuities. This sum was relatively low for the significant ones. This was caused by the vertical profile of discontinuity occurrence in the troposphere. From 300 hPa down to 500 hPa (Figure 21), the sum was high, but there were regions where this sum was much lower (over Africa and the Indian Ocean, over China, and over South America and the South Atlantic).

**Figure 15.** Geographical distribution of all (**upper panel**) and the significant (**lower panel**) discontinuities over the globe at 0.3 hPa.

**Figure 16.** The same as Figure 15, but for 0.4 hPa.

**Figure 17.** Geographical distribution of all (**upper panel**) and the significant (**lower panel**) discontinuities over the globe at 1 hPa.

**Figure 18.** The same as Figure 17, but for 5 hPa.

**Figure 19.** Geographical distribution of all (**upper panel**) and the significant (**lower panel**) discontinuities over the globe at 30 hPa.

**Figure 20.** The same as Figure 19, but for 300 hPa.

**Figure 21.** Geographical distribution of all (**upper panel**) and the significant (**lower panel**) discontinuities over the globe at 500 hPa.

#### **4. Discussion**

When we look at the ozone concentration patterns at the uppermost model layers, we see the strange patterns down to 0.5 hPa (Figures 1 and 2). According to our opinion, this is caused by the model itself, and these patterns are not real. So these layers are not suitable for trend analyses. The other groups of patterns are the ones with maximal concentration over the subpolar latitudes (from 0.7 to 2 hPa) of the winter hemisphere. It is very hard to explain this behaviour, because one can expect the main influence of solar activity to ozone at these heights, and in winter, there is very little solar radiation at the latitudes of maximum occurrence. In addition, at these heights, the ozone concentration is under the control of photochemistry, and the Brewer–Dobson circulation cannot act on them strongly. So one possible explanation is again the internal model variability, and the exploitation of these layers in trend analyses is problematic. Below these heights, the patterns with equatorial maximum and polar minimum were observed (from 3 down to 10 hPa). This is caused by the direct influence of solar radiation on ozone formation, which maximizes ozone in low latitude, and the Brewer–Dobson circulation is weak at these levels [24]. This upper stratospheric pattern can be explained by the theory, so in these heights, MERRA 2 gives reasonable results. From 30 hPa down to the tropopause, the maximal concentration was seen over high latitudes of both hemispheres and the minimal one over the equator. This is caused by the Brewer–Dobson circulation, which transports ozone rich air to high latitudes [25,26]. It is interesting these patterns persist deep into the troposphere, especially in the Northern Hemisphere, which may be caused by stratosphere–troposphere exchange in the vicinity of the Aleutian low.

The other topic of this paper, which is worth discussing, is the temporal occurrence of discontinuities in the period 1980–2017. In the uppermost model layer, the maximum of their occurrence was seen in 1993. One of the explanations of 1993 maximum is the Pinatubo eruption in 1991, but this maximum occurred only in the uppermost model layers. [27] found the enhanced wave activity in the two winters just after the Pinatubo eruption in similar heights. Thus the effect of the Pinatubo eruption may be real in the uppermost model layers. Going down this maximum disappeared, and the new one set up in 2004, and it was strongest at 1 hPa. This maximum can be clearly explained by the transition from SBUV to EOS Aura data in 2004 [12,28]. Below this layer, the flat temporal distribution of discontinuities was seen down to the troposphere, so no maximum in 2004 was observed. We can explain this fact by the following reasons. We can speculate the Pettitt test was either not suitable for this analysis, or the discontinuities in 2004 were not the most important in time series. The Pettitt test is widely used in atmospheric research, so we have no reason why this test does not work well for MERRA 2 ozone concentration data. If the 2004 discontinuities were not so huge, this fact could open the possibility of using MERRA 2 ozone data in trend analysis. The decision regarding these two hypotheses will be a matter of future research. In the troposphere, we observed wide maximum again in 1993, but not in 2004. In future work, we shall look at the wave activity in the upper troposphere and search its connection with the discontinuities occurrence peak in 1993.

Shangguan et al. [12] also detected discontinuities in 2015 when version 4.2 of the Microwave Limb Sounder (MLS) ozone data were used instead of version 2.2 and in 1998 when the Advanced TIROS (The Television Infrared Observation Satellite) Operational Vertical Sounder (ATOVS) was launched. In our paper, these discontinuities were not detected due to these reasons: We used data in the period of 1980–2017, so the year 2015 is at the edge of the series, and the Pettitt test could not capture this year. In each grid, we took into account only one main discontinuity, so the year 1998 could not be the main discontinuity, and thus, we did not notice it.

When we look at the vertical profile of discontinuity occurrence, we see the two main areas: the uppermost model layers and the tropospheric layers. From about 2 hPa down to the troposphere, we observed fewer discontinuities. High occurrence of discontinuities in the lowermost stratosphere and the upper troposphere can be explained by the lower quality of SBUV ozone data in these layers [29]. The share of the significant discontinuities was high in the upper model layer, while this share was low in the troposphere. The less number of (significant) discontinuities, the better for trend analysis.

We can conclude that due to ozone concentration patterns and the occurrence of discontinuities, the layers above 4 hPa are not suitable for trend analyses in MERRA 2. Below 4 hPa, the MERRA 2 gives explainable results expected from the theory, which also promotes using these data in the trend analysis. According to our results, the transition from SBUV to EOS Aura data does not have a major impact on discontinuity occurrence at these layers. The worst condition for trend analyses was the occurrence of the significant discontinuities in time series, and the share of grids with them was relatively low below 4 hPa. The better situation was the existence of discontinuities, which were not significant. Their number was higher than the significant ones. The grids with no discontinuities were not rare in the MERRA 2 ozone concentration data.

We must also deal with the possibility of the occurrence of discontinuities with natural origin. They must be included in trend analysis because they are real. It is very difficult to distinguish between artificial and real discontinuities because the real ones occur when there is some natural change in the atmosphere. We must consider the possibility of atmospheric changes connected to the change in ozone content and with increasing concentration of CO2. Sudden water vapour drops in the stratosphere around 2000. In ref [30] can also be the manifestation of the changes in the atmospheric behaviour.

Almost all the modern reanalyses had a temporal discontinuity in 1998 when the ATOVS observations became available, and the reanalyses either switched immediately or transitioned from the TIROS Operational Vertical Sounder (TOVS) to the ATOVS over several years. These discontinuities were not consistent in each latitude or pressure level. The disagreement in the polar latitudes of the southern Hemisphere extended lower into the stratosphere than in the northern hemispheric ones. The zonal wind variance was smaller than the temperature variance in the polar latitudes, but had a similar temporal difference between the TOVS and ATOVS time periods. For the tropics, the zonal wind variance is much larger than in the polar regions, and the disagreement between the semi-annual oscillation (SAO) and QBO zonal winds is quite large. Thus, improving equatorial winds in future reanalyses is an important goal. MERRA 2 reanalysis also showed discontinuities in 2004 when MLS observations were included. The disagreement was evident mainly for temperature or zonal wind. More details about discontinuities in reanalyses can be found in ref [31].

#### **5. Conclusions**

The main achievements concerning the ability to utilize the MERRA 2 ozone data for trend studies based on the Pettitt test results and on the analyses of monthly averages are the following:


**Author Contributions:** This paper was done in close collaboration of the authors, but substantial part of the works was done by P.K.

**Funding:** Support by the Grant Agency of the Czech Republic via Grant 18-01625S is acknowledged.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*
