Technical note: a simple procedure for removing temporal discontinuities in ERA-Interim upper stratospheric temperatures for use in nudged chemistry-climate model simulations

. This note describes a simple procedure for removing unphysical temporal discontinuities in ERA-Interim upper stratospheric global mean temperatures in March 1985 and August 1998 that have arisen due to changes in satellite radiance data used in the assimilation. The derived temperature adjustments (offsets) are suitable for use in stratosphere-resolving chemistry-climate models that are nudged (re-laxed) to ERA-Interim winds and temperatures. Simulations using a nudged version of the Canadian Middle Atmosphere Model (CMAM) show that the inclusion of the temperature adjustments produces temperature time series that are devoid of the large jumps in 1985 and 1998. Due to its strong temperature dependence, the simulated upper stratospheric ozone is also shown to vary smoothly in time, unlike in a nudged simulation without the adjustments where abrupt changes in ozone occur at the times of the temperature jumps. While the adjustments to the ERA-Interim temperatures remove signiﬁcant artefacts in the nudged CMAM simulation, spurious transient effects that arise due to water vapour and


Introduction
Stratosphere-resolving chemistry-climate models (CCMs) are powerful tools for simulating and understanding the impacts of ozone depletion/recovery and increasing greenhouse gas emissions on the climate system and for predicting its future changes. However, because free-running CCMs are unable to reproduce the observed meteorology, the simulated fields cannot be compared to observations on a day-to-day basis. Furthermore, biases in one model variable can have knock-on effects on other variables. For example, a bias in temperature in the upper stratosphere (due to a missing process in a model's radiative transfer scheme) can impact on chemistry, thus causing a bias in another variable such as ozone, even if the chemistry scheme were perfect.
Recently, there has been a concerted effort to nudge CCMs to the horizontal winds and temperatures from reanalyses using Newtonian relaxation. Constraining the dynamical fields to follow the observations through nudging enables day-today comparisons of other model variables to observations to be made and biases in the model "physics" and chemistry to be identified. To this end the International Global Atmospheric Chemistry (IGAC) and Stratospheric Processes And their Role in Climate (SPARC) Chemistry-Climate Model Initiative (CCMI) has proposed a set of specified-dynamics simulations in which the CCMs are nudged to reanalysis data sets from 1980 to 2010 (Eyring et al., 2013). The proposed REF-C1SD nudged simulations will be used for upcoming ozone and climate assessments.
This approach inherently relies on the reanalysis used to do the nudging not containing spurious temporal inhomogeneities. This is a challenge for reanalyses because of the inhomogeneous nature of the observations. The reanalysis centres have made major advances in dealing with these issues through bias correction, with the European Centre for Medium-Range Weather Forecasts (ECMWF) "Interim" Reanalysis (ERA-Interim; Dee et al., 2011) being the first reanalysis to use fully automated variational bias corrections of satellite observations. However, bias correction relies on having multiple data sets, in particular "anchoring" data sets such as radiosondes. Unfortunately, in the upper stratosphere, where there is only one source of long-term temperature data (operational nadir-viewing satellites), bias correction cannot be effectively applied if, as in the case of ERA-Interim, model bias is large and adjustment towards a biased model state is to be avoided. Dee and Uppala (2008) discuss the difficulties associated with assimilating radiances from the Stratospheric Sounding Unit (SSU) and the Advanced Microwave Sounding Unit (AMSU-A). They show that the transition from SSU to AMSU-A data in August 1998 introduces large global mean temperature increments in the upper stratosphere in ERA-Interim (their Fig. 18), which result in a sudden jump in global mean temperature in the reanalysis, which is largest at 1 hPa (their Fig. 20). If that jump is removed, a slight cooling trend in the upper stratosphere is observed over the period they analysed from 1989 to 2005.
Here we discuss a simple procedure for removing unphysical temporal discontinuities in global mean ERA-Interim temperature data in the upper stratosphere for use in nudged CCMs. Corresponding adjustments to the reanalysis wind fields are not required since to a very good approximation changes in global mean temperature do not impact on dynamics in the stably stratified stratosphere. The temperature adjustment proposed here will be useful for high-top CCMs used for the REF-C1SD simulations for CCMI. We apply this procedure to the jump in August 1998 discussed in Dee and Uppala (2008) and also to another large jump in March 1985, which was outside of their analysis period. This 1985 jump is due to a change in bias of SSU-3 data from NOAA-9 and subsequent satellites, for which pressure in the instrument's carbon dioxide cell was higher than on earlier satellites (Kobayashi et al., 2009).
We provide a set of equations and a tabulated set of coefficients that can be easily used to generate the global mean temperature adjustments. We also provide a link for obtaining data files containing those adjustments. We must stress that our adjustments do not correct the ERA-Interim temperatures, they only remove the unphysical temporal discontinuities, which are an artefact of the assimilation process.
We assess the impact of including the temperature adjustments in a nudged version of the Canadian Middle Atmosphere Model (CMAM). The simulation without the temperature adjustments exhibits large discontinuous changes in global mean temperature in the upper stratosphere, which are absent when the adjustments are applied. The inclusion of the adjustments is shown to have a large impact on the simulated upper stratospheric ozone, which exhibits large jumps contemporaneous with the temperature jumps when the unadjusted nudging temperatures are used.

Adjusting the ERA-Interim global mean temperatures
This section describes the methodology used for removing the temporal discontinuities (jumps) in the ERA-Interim temperatures in the upper stratosphere. We consider only the global mean since, to a first approximation, it is unaffected by dynamics and is therefore close to radiative equilibrium (e.g. Fomichev, 2009). Since global mean temperature varies slowly in time (due to changes in radiative forcing), the jumps can be identified in an unambiguous manner. We start by identifying the times at which the jumps occur and then describe the procedure for computing and applying the adjustments to the data, which we shall refer to as offsets. The offsets are computed from monthly mean data, but are interpolated to the 6-hourly data that are used to nudge the model. Figure 1 shows global and monthly mean temperature time series for the top five pressure levels of ERA-Interim. Anomalies computed with respect to the long-term  monthly means are shown here to highlight the discontinuities. Superimposed on the overall cooling trend at all levels are large temperature jumps in March 1985 and August 1998, denoted by the thin vertical lines. These are the two jumps that will be removed. The jumps occur simultaneously at all levels, and decrease in magnitude with increasing pressure. The jump in August 1998 results from the switch from SSU to AMSU-A satellite radiances (Dee and Uppala, 2008;Kobayashi et al., 2009). It is not instantaneous, but is spread out over a period of about 10 days, as seen by the red curves in Fig. 3b, which shows the 6-hourly data at 1 hPa. The jump in March 1985 exhibits a much different temporal behaviour, for not only is there a shift in the mean value, but there is also a spike just before the jump, as is more clearly seen in Fig. 3a (red curves). This complex temporal behaviour is caused by varying SSU-3 data counts from the NOAA-7 and NOAA-9 instruments, which as noted earlier have different biases, during their several-month overlap period (A. Simmons, personal communication, 2013). There is also a jump in 1979 associated with the start of assimilation of SSU-3 data, which we do not attempt to adjust since the ERA-Interim record prior to this time is very short. Figure 1 shows that at 7 hPa the jumps in 1985 and 1998 are not readily distinguished from the geophysical noise. Consequently we shall adjust only the top four pressure levels, namely 1, 2, 3, and 5 hPa. We note that our procedure does not remove possible slow drifts in the reanalysis data that might arise from long-term changes in the satellite data. Figure 2 illustrates the steps used in computing the monthly mean offsets. We focus on 1 hPa since that is where the temperature jumps are largest; the methodology is identical for the other three pressure levels. Figure 2a shows the original ERA-Interim global mean temperature time series, which differs from the top curve in Fig. 1 in that it includes the annual cycle. These data are deseasonalized by subtracting off the fit to the first four sine and cosine harmonics of the annual cycle. The fits are calculated independently for three time periods, namely 1980 to 1984, 1986 to 1997, and 1999 to 2011, which exclude the years when the jumps occur. The resulting deseasonalized time series is given by the black curve in Fig. 2b. Since the time axis starts in 1980, the jump in 1979 is not seen here. The deseasonalized time series is then fit to a function that includes a linear trend, the monthly mean F10.7 solar flux time series recorded by the Geological Survey of Canada (http://www.spaceweather. ca/data-donnee/sol_flux/sx-5-eng.php) and the global mean stratospheric aerosol optical depth time series updated from Sato et al. (1993) as a proxy for volcanic activity (http:// data.giss.nasa.gov/modelforce/strataer/). The linear trend accounts for cooling due to the increasing concentrations of CO 2 (and, before 1998, ozone loss), while the solar flux term accounts for solar cycle effects, which have a significant impact on upper stratospheric temperatures. The volcanic term is included because it improves the fit, although the upper stratosphere is far from the main region of volcanic influence; omitting the stratospheric aerosol optical depth term has only a minor effect on the values of the offsets. (In response to a comment by one of the reviewers, we have also tested fitting a term to account for the quasi-biennial oscillation but found that there was little improvement in the fit nor effect on the value of the offsets, as expected for global mean temperatures.) The fits are computed separately over the three periods given above, and the time series of the fits are then extended across the years in which the jumps occur. For example, data from January 1980 to December 1984 were used to construct the fit for the first period; the fit was then extended forwards in time to cover 1985. Similarly, the second period was fit using data from January 1986 to December 1997, and the fit was extended backwards in time to cover 1985. The same procedure was used to have the second and third periods span 1998. The red curves in Fig. 2b show the fitted data, extended slightly beyond the ends of the data used in computing the fits in order to show the overlap region that is used to calculate the annual mean offsets, which are evaluated at the mid-point of the jump year denoted by the red crosses.
A simpler approach for generating the offsets would be to compute global mean temperature differences for specific periods before and after the jumps. The problem with that approach, however, is that if the time periods are too long, the cooling trend aliases into the average; while if the periods are too short, the systematic bias is harder to characterize. Our approach mitigates both of these problems and produces annual mean offsets that yield a smoother transition over the jumps.  At this stage of the analysis, a decision has to be made as to how to apply the offsets. Since there are no independent observations that can be used to validate the global mean temperatures, we will apply the offsets to only the first and second time periods, leaving the post-1998 period as the (unaltered) reference, noting that referencing the adjustment procedure to the last period is an arbitrary choice. Figure 2c shows the deseasonalized time series with the annual mean offsets applied (only the 1998 offset to the middle segment of the data set, and the sum of the two offsets to the first segment of the data set), which with the inclusion of the annual cycle yields Fig. 2d. (The precise way that these offsets are applied is discussed below.) The final step involves the inclusion of the monthly mean offsets. As the post-1998 ERA-Interim data were defined as the absolute reference for the record, seasonal cycle offsets were calculated as the difference between the fitted seasonal cycle of the ERA-Interim monthly temperatures for a particular period and the fitted seasonal cycle for the post-1998 period. The monthly mean offsets for a particular period are then the sum of the seasonal cycle offsets, which by construction have an annual mean of zero, and the annual mean offset calculated from the fitting described above. The final time series with the monthly mean offsets applied is shown in Fig. 2e. Note that the spike in March 1985 remains since there is no way to remove it using our procedure.
The end product of this part of the analysis procedure is two sets of monthly mean offsets for each of the four pressure levels. It now remains to interpolate those offsets to a 6 h grid and to smoothly connect them over the month where the jumps occur. The former is achieved by fitting the monthly where D = 365/(2π ), d the day of the year (e.g. 0, 0.25, 0.5, 364.75 for 6 h data, with leap years excluded), l the pressure level (1, 2, 3, or 5 hPa), n the time period (1 or 2), and m the harmonic of the annual cycle. The values of the fit coefficients A m l,n , B m l,n , and C l,n are given in Table 1. The final step is to smoothly join the 6-hourly offsets for the adjacent periods, because the jump between periods is not instantaneous (see Fig. 3). This is achieved by applying an exponential function of the form E l,n (d) = exp (−(d − d n )/τ l,n ), where d n is the day of the year of the jump for each period (e.g. d 1 = 80.5 for 12:00 UTC 21 March), and τ l,n is an e-folding timescale. The latter is specified as τ l,1 = 3 days for all levels, τ l,2 = 5 days for l = 1 hPa, and τ l,2 = 2 days for the remaining levels. These values were chosen to yield a smooth transition, judged by eye.
Denoting the offsets at level l as T l , we have for t < t 1 (where t is time and t 1 denotes the first "jump time", 12:00 UTC 21 March 1985) and d is, as before, the day of year: T l (t) = T * l,1 (d). ( In order to smoothly merge the exponential function that is applied immediately after the first jump with the 6-hourly offsets for the second period, we use for t 1 ≤ t < t 1 + 25 days: Atmos. Chem. Phys., 14, 1547-1555, 2014 www.atmos-chem-phys.net/14/1547/2014/ Table 1. Coefficients A m l,n , B m l,n and C l,n used in Eq. (1) for generating the 6-hourly global mean temperature offsets using the first four harmonics of the annual cycle for the two time periods. Period n =1 extends up to 12:00 UTC 21 March 1985; period n = 2 extends from 12:00 UTC 21 March 1985 to 00:00U TC 1 August 1998. Units are K. In order to provide more precision, the coefficients A m l,n and B m l,n are multiplied by 10.

Period (n) Pressure (l)
C l,n A m l,n × 10 B m l,n × 10 For t 1 + 25 ≤ t < t 2 (where t 2 denotes the second jump time, 00:00 UTC August 1, 1998), and for t 2 ≤ t < t 2 + 25 days, Finally, for t ≥ t 2 + 25 days, Figure 3, which shows a blow-up of the jumps in 1985 and 1998, helps to illustrate the use of the above set of equations. The red curves denote the original 6 h ERA-Interim data at 1 hPa; the blue curves are the adjusted data (i.e. with the offsets T l added on). The jump times t 1 and t 2 are denoted by the left-hand thin vertical lines. Figure 4a shows the offsets T l at all four pressure levels computed using Eqs. (1)-(6) and monthly averaged for plotting purposes. As expected, the largest component of the calculated offset is the annual mean. The corresponding global and monthly mean temperature anomalies with and without the offsets included (the adjusted and unadjusted anomalies, respectively) are given by the black and red curves, respectively, in Fig. 4b. The red curves are identical to those shown in Fig. 1.

Model description and simulations
The Canadian Middle Atmosphere Model is a chemistryclimate model that extends up to the lower thermosphere (∼ 100 km). The version used here is nudged to the 6-hourly horizontal winds and temperatures from ERA-Interim. The nudging is performed in spectral space and is applied only to horizontal scales with n ≤ n max = 21, where n is the total wave number. The nudging tendency has the form − (X − X R )/τ r , where τ r = 24 h is the relaxation timescale and X and X R are, respectively, the model and reanalysis nudged prognostic variables, namely vorticity, divergence and temperature (on model levels). Above 1 hPa the model evolves freely. The nudging configuration is identical to that described in McLandress et al. (2013), with the exception that here we use a constant value for τ r at all heights. Additional information about the free-running version of CMAM can be found in Scinocca et al. (2008). Three simulations are performed. The first is a nudged simulation that includes the global mean offsets, T l , to the ERA-Interim temperatures computed using Eqs. (1)-(6). (Since the model is spectral, the offsets are added to the global mean harmonic. In a grid-point model the same offset would be applied at each latitude and longitude.) The second is identical to the first but without the offsets (i.e. using the original ERA-Interim temperatures). These two simulations will be referred to as the "adjusted" and "unadjusted" runs. The third is a free-running simulation (i.e. without any nudging), which is otherwise identical to the two nudged runs; it will be referred to as the "free run". The simulations extend from 1979 to 2010, and are initialized on 1 January 1979 using an earlier run nudged to the ERA40 (Uppala et al., 2005) winds and temperatures. Monthly and annually varying sea surface temperatures and sea-ice distributions are prescribed using observations (Rayner et al., 2003). The radiative forcings and chemical boundary conditions are the same as those used in SPARC CCMVal (2010). Solar variability is not included in any of these runs.

Model results
Here we examine the impact of the ERA-Interim global mean temperature adjustments on the nudged model results.
We focus on the global mean response since it is the firstorder effect of including the adjustments. We also discuss the simulated ozone field since it is strongly tied to temperature in the upper stratosphere. Since the ozone is freely evolving (i.e. unnudged), differences between the ozone for the adjusted and unadjusted runs demonstrate how the temperature in the nudged region of the model affects chemical constituents. Results from the free run serve to highlight the model biases. Figure 5a shows the deseasonalized global and monthly mean temperatures at 1 hPa for the three simulations shown here by the solid curves (the dotted curve will be discussed shortly). As expected the results for the adjusted and unadjusted runs closely resemble the ERA-Interim time series shown in Fig. 2b and c, with the unadjusted run (red) exhibiting the jumps in 1985 and 1998 and the adjusted run (blue) having a continuous cooling trend void of the two jumps. The free run (black solid) has a warm bias of at least 3 K relative to ERA-Interim, which is presumably related to a bias in the CMAM radiation code (see Chapter 3 of SPARC CCMVal, 2010). The corresponding deseasonalized global mean ozone is shown in Fig. 5b. Since ozone is under photochemical control at 1 hPa, the differences between the three solid curves are primarily attributed to differences in temperature. Prior to 1998 the unadjusted run has the highest ozone concentration, since the lower temperatures in this time period result in reduced photochemical destruction of ozone and hence more ozone. The sudden temperature increase in 1998 in that run produces an immediate decrease of ozone. The inclusion of the temperature adjustments results in a smooth transition of ozone across 1998. The impact of the 1985 temperature jump on ozone is more difficult to see since ozone happens to peak about that time.
The temperature and ozone trends for the free run (solid black curves in Fig. 5) are of opposite sign before and after 1985, which at first seems odd since in the early 1980s one would expect to see cooling in response to the increasing concentrations of CO 2 and ozone loss in response to the increasing concentrations of ozone-depleting substances. That this does not occur is due to the initial conditions, which as stated earlier are from a run nudged to ERA40 data. The run nudged to ERA40 had significantly higher concentrations of stratospheric water vapour in 1980 than in a corresponding free-running version of the model (not shown) because of higher tropical tropopause temperatures. Therefore, the transient adjustment of the water vapour from its high initial values to the lower values that are characteristic of the free-running model results in a warming (i.e. less radiative cooling; see Fomichev, 2009), which is strong enough to counteract the global mean cooling due to increasing concentrations of CO 2 . Note that it takes about 2 yr for this signal to propagate from the tropical lower stratosphere up to the stratopause, through the Brewer-Dobson circulation; see ted black) results from a slightly older version of the freerunning model that was initialized using an unnudged run (i.e. a run that did not have this moist "bias"). In this run (denoted here as the "old free run" to distinguish it from the "free run"), cooling commences from the start of the plotted time series, and by 1985 the temperature closely tracks that of the free run (solid black), with both exhibiting cooling in response to the increasing concentration of CO 2 . The initial conditions also have an important impact on global mean ozone at 1 hPa, as seen by comparing the results for the two free runs in Fig. 5b. The ozone increase in the early 1980s in the free run is again a result of the moist initial conditions: as the model adjusts to the drier conditions inherent in the free-running model, the photochemical destruction of ozone through the HO x cycle is reduced, resulting in increasing concentrations of ozone. (The effect on ozone cannot be driven primarily by temperature, which would act in the opposite direction.) Figure 4b of Shepherd (2002) shows the roughly 2-year delay in the onset of the upper stratospheric HO x perturbation following a change in tropical tropopause temperatures. By 1985 this transient adjustment abates, and ozone starts to decline in response to increasing concentrations of ozone-depleting substances (ODSs), decreasing up until the late 1990s, after which it starts to increase due to decreasing concentrations of ODSs.
This transient dependence of upper stratospheric ozone on the initial conditions is also seen in the two nudged runs (Fig. 5b), which both show increases from 1980 to 1985, much like in the "free run" results. Here the spurious transient arises from the difference in stratospheric water vapour in the run nudged to ERA40, which provide the Atmos. Chem. Phys., 14, 1547-1555 initial conditions, and the runs nudged to ERA-Interim, with the latter being significantly drier (results not shown). A similar transient adjustment can be expected in any nudged simulation, and thus the upper stratospheric ozone behaviour prior to 1985 in the CCMI REF-C1SD simulations should probably be disregarded unless the nudged climatology is identical to the spin-up climatology.
The vertical extent of the impact of the temperature adjustments is shown in Fig. 6, which shows profiles of global mean temperature averaged over the period between the temperature jumps in 1985 and 1998. Not surprisingly, the largest impact is at 1 hPa, where the adjusted run is about 2.5 K warmer than the unadjusted run. The difference changes sign at about 3 hPa, with the adjusted run being about 1 K colder at 4 hPa, in accordance with the change in sign of the temperature adjustments at these levels (Fig. 4a). With the exception of the 1 hPa level for the adjusted run, the nudged runs are slightly warmer than their ERA-Interim counterparts (dashed and dotted curves). This is due to the fact that at all heights plotted here the free run is warmer than ERA-Interim, and the nudging pulls the model temperatures from their free-running values towards the reanalysis.
The final two figures show the full seasonal cycle of temperature and ozone at 1 hPa (Fig. 7) and 2 hPa (Fig. 8). At 2 hPa the impact of including the adjustments is smaller than at 1 hPa, as expected based on Fig. 4a.

Summary
This technical note has described a simple procedure for removing unphysical temporal discontinuities in ERA-Interim global mean temperatures in the upper stratosphere in March 1985 and August 1998 for the purpose of improving chemistry-climate model simulations that are nudged to ERA-Interim winds and temperatures. The temperature adjustments are derived from the global and monthly means at 1, 2, 3 and 5 hPa and are interpolated to a 6 h grid. Adjustments are not made to the post August 1998 data since they are treated as the absolute reference. We also note that our procedure cannot identify slow drifts that might exist within each of our analysis periods.
Two simulations using a version of the Canadian Middle Atmosphere Model that is nudged to the three-dimensional 6 h ERA-Interim horizontal winds and temperatures are www.atmos-chem-phys.net/14/1547/2014/ Atmos. Chem. Phys., 14, 1547-1555, 2014 performed in order to assess the impact of including the temperature adjustments. The global mean temperatures in the upper stratosphere in the simulation without the adjustments exhibit the large temperature jumps found in ERA-Interim. Moreover, the ozone time series in this simulation also exhibit large discontinuities at the times of the temperature jumps. This occurs because ozone is under photochemical control in the upper stratosphere. The inclusion of the adjusted nudging temperatures produces global mean temperatures and ozone that are devoid of the spurious jumps. We encourage other modellers who are planning on performing nudged chemistry-climate model simulations to use the adjustment procedure described here. To this end we have also made available both ASCII and netCDF files containing 6 h time series of the global mean ERA-Interim temperature adjustments (offsets) from 1979 to 2010 at 1, 2, 3, and 5 hPa. These data can be found at http://www.cccma.ec.gc.ca/data/ cmam/cmam30/era_interim_adjustment.
An important side issue arising from this work that is relevant to the REF-C1SD simulations for CCMI is the initialization of simulations nudged to ERA-Interim from 1979 onwards. We have initialized ours using a spin-up nudged to ERA40. However, since the tropical tropopause in ERA40 is warmer than in ERA-Interim, stratospheric water vapour in runs nudged to ERA40 is systematically higher than in runs nudged to ERA-Interim. This leads to a moist bias in the initial condition for the simulation nudged to ERA-Interim, which takes about 5 yr to disappear, consistent with the mean age of air in the upper stratosphere. During this adjustment period of decreasing water vapour, there is an associated decrease in the destruction of ozone through HO x chemistry, which outweighs ozone loss due to increasing concentrations of ozone-depleting substances and leads to an otherwise puzzling (and physically spurious) increase in upper stratospheric ozone during the first 5 yr of the ERA-Interim period. It is only after 1985 that upper stratospheric ozone starts to decline in the nudged simulation. Thus, chemical constituents with a strong dependence on stratospheric water vapour (e.g. upper stratospheric ozone) should probably be disregarded from 1979 to about 1985 in simulations nudged to ERA-Interim unless the nudged climatology is identical to the spin-up climatology.