Long-term Midlatitude Mesopause Region Temperature Trend Deduced from Quarter Century (1990–2014) Na Lidar Observations

The long-term midlatitude temperature trend between 85 and 105 km is deduced from 25 years (March 1990– December 2014) of Na Lidar observations. With a strong warming episode in the 1990s, the time series was least-square fitted to an 11-parameter nonlinear function. This yields a cooling trend starting from an insignificant value of 0.64 ± 0.99 K decade −1 at 85 km, increasing to a maximum of 2.8 ± 0.58 K decade −1 between 91 and 93 km, and then decreasing to a warming trend above 103 km. The geographic altitude dependence of the trend is in general agreement with model predictions. To shed light on the nature of the warming episode, we show that the recently reported prolonged global surface temperature cooling after the Mt Pinatubo eruption can also be very well represented by the same response function .

1 Introduction Roble and Dickinson (1989) estimated the effects of hypothetical future increases in greenhouse gas concentrations on the global mean structure and predicted considerable cooling in the mesosphere and thermosphere.About this time, a number of long-term temperature observations in the mesopause region (80-110 km) were initiated or reinitiated at locations in the Northern Hemisphere with passive OH emissions and/or active probes, such as Na lidar and falling spheres.These observations and those in the Southern Hemisphere via OH emission, as well as the long series of Russian rocket measurements and OH emissions between about 1960 and 1995 over a wide range of latitudes, measured cooling trends in the mesopause region ranging from 0 to ∼ 10 K decade −1 , suggesting that after 2 decades, the observed trend remains uncertain (Beig, 2006).These observational temperature trend results were referenced in Table I of She et al. (2009).
Based on the nocturnal lidar temperatures acquired between March 1990 and December 2007 (data set (90-07)), the same paper reported a linear long-term trend, starting from an insignificant cooling trend of 0.28 ± 1.32 K decade −1 at 87 km, reaching a maximum value of 1.55 ± 1.15 K decade −1 at 91 km, and turning into a warming trend above 102 km.The magnitude and altitude dependences are consistent with the prediction of the Spectral Mesosphere/Lower Thermosphere Model (SMLTM) (Akmaev et al., 2006) and of the Hamburg Model of the Neutral and Ionized Atmosphere (HAMMONIA) (Schmidt et al., 2006).Subsequent substantial reviews on thermospheric trends, Lašttovička et al. (2012) and Cnossen (2012), included some studies on the mesopause region neutral temperatures.Recent observational reports on mesopause region temperature trends at specific altitudes include Offermann et al. (2010) and Hall et al. (2012), based on ∼ 10year data sets.The former utilized the annual mean OH imager temperatures between 1988 and 2008 over Wippertal (51 • N, 7 • E) and reported a long-term trend at 87 km of −2.3 ± 0.6 K decade −1 in temperatures; the latter utilized meteor wind radar between October 2001 and October 2012 at Svalbard (78 • N, 16 • E) calibrated by satellite measurements and reported a temperature trend at 90 km of Published by Copernicus Publications on behalf of the European Geosciences Union.

Na Lidar data sets and long-term regression analysis
The Colorado State University (CSU) Na lidar performed nocturnal mesopause region temperature observations between March 1990 and March 2010 at Fort Collins, CO (41 • N, 105 • W).It employed a vertical beam between 1990 and 2001.Since 2002, the lidar has operated in 2-or 3-beam geometry for simultaneous temperature and wind measurements, leading to two or three mean temperatures at a given altitude each night.The lidar was relocated to Utah State University (USU) and has continued its regular observation at Logan, Utah (42 • N, 112 • W) since September 2010.Because of similar geographical coordinates, we combine the data from both locations to form a data set from March 1990 to December 2014, denoted as (90-14_Avg).The extension _Avg indicates that, unlike the data set employed in previous publications such as (90-07), which utilized temperatures from the beam with the largest signal at 3.7 km vertical resolution, here we use the average of temperatures acquired from two or three beams at 2.0 km vertical resolution.
As an overview, we plot the 25 years of nightly mean temperatures at 86 km, which shows large annual and semiannual variation, and at 99 km, an altitude with minimal annual and small semiannual variation (She and von Zahn, 1998), respectively, in Fig. 1a and b.The data acquired at CSU (March 1990-March 2010) is in black and that acquired at USU (September 2010-December 2014) is in blue; apart from a small data gap in 2010, the two sets of data blend nicely.From Fig. 1a summer is about 60-80 K cooler than winter at 86 km.At 99 km one can see long-term temperature variation.The 81-day averaged daily F10.7 solar flux also plotted in the figure, in the red curve, shows that the nightly mean temperatures track the variation in solar flux after 1993.Note that there exists a warming episode after the Mt Pinatubo Eruption (MPE), t MPE = 1.45 years, which peaked in 1993 and mostly died away near the beginning of 1999 for altitudes between 88 and 102 km (Fig. 3a).Since the warming episode is in our data, we must account for it in the analysis, whether its causes are fully understood or not.As a result, a nonlinear least-square regression analysis is required for long-term study.
Following She et al. (2009), who performed regression analysis on a shorter data set (90-07), a time series with 894 points, we express the nocturnal temperature at each altitude, T (z, t), of (90-14_Avg), a time series of 1200 points, as where, T fit (z, t) = α(z) + A 1 (z) cos(2π t) where t is time in years from 1 January 1990, α(z) is independent of time, and the 4 A-B terms represent annual and semiannual variations.The three long-term effects have the amplitudes β(z), γ (z), and δ(z), in this model.The strong warming episode in our data, initially attributed to the Mt Pinatubo eruption in June 1991 (She et al., 1998), is represented by an amplitude γ (z) times P (z, t) = 2/ {exp(t 0 − t)/t 1 + exp(t − t 0 )/t 2 }, with parameters t 0 (z), t 1 (z), and t 2 (z), for the delay, rise, and decay time, respectively.The delay time here is relative to 1 January 1990, with the Mt Pinatubo eruption at 1.45 years (see Fig. 1b).Other long-term responses include δ(z), the solar response in K/SFU with Q 81 (t) being the 81-day averaged F10.7 solar flux, and β(z), the linear trend in K years −1 .The residual from the best fit is T Res (z, t).Since all effects of comparable strengths must be included in the time series for the nonlinear regression analysis (Akmaev, 2012) and the warming episode, solar activity and linear trend are not independent, the best fit of one term will affect that of the other and they will depend upon the length of the data set.

Temperature trend deduced from quarter century lidar data
The long-term linear trend of the 11-parameter fit to the long data set, F-11P(90-14_Avg), is shown in Fig. 2 along with F-7P(90-14_Avg), deduced from the 7-parameter fit by setting γ (z) = 0. Also shown are the published results from F-11P(90-07) and F-7P(90-07) from the shorter data set from March 1990 to 2007.As expected, the uncertainty from the 25-year data set is smaller than that from the 18-year data set.The cooling trend in F-11P(90-14_Avg) starts from an insignificant value of 0.64 ± 0.99 K decade −1 at 85 km, increases to a maximum of 2.8 ± 0.58 K decade −1 between 91 and 93 km, and then gradually decreases to a warming trend above 103 km.Turning from a cooling to a warming trend above 100 to 120 km in geographic altitudes is the result of the cooling and contraction of underlying atmosphere, the lower thermosphere, the mesosphere, and the stratosphere; it is predicted by models (Akmaev, 2012;Qian et al., 2013).This does not occur in the seven-parameter analysis (see F-7P(90-14_Avg) in Fig. 2).A similar difference between F-11P(90-07) and F-7P(90-07) is also evident.To our knowledge, metal resonance lidar is the only ground-based instrument that covers the entire mesopause region, and ours is the only Na lidar with a long enough data set to see this trend reversal.
Compared to the trends deduced from the shorter data set, (90-07), we note that the difference between F-7P(90-07) and F-11P(90-07) is bigger than the difference between F-7P(90-14_Avg) and F-11P(90-14_Avg) because the influence of the warming episode on the temperature trend is reduced in a longer data set.Statistically, the results from the longer data set are more accurate; the mean uncertainty between 88 and 102 km is 0.6 and 1.3 K decade −1 , respectively, for F-11P(90-14_Avg) and the previously published F-11P(90-07).However, below we investigate the discrepancy between the two 11-parameter analyses, i.e., the 25-year data set has a larger cooling trend by ∼ 1 K decade −1 .
Since the three long-term effects with magnitudes β(z), γ (z), and δ(z) are not independent in our analysis, we can understand their mutual influences by realizing that, in addition to the annual and semiannual variations, the observed temperature at a given time is the sum of three contributions β(z)t, γ (z)P (z, t), and δ(z)Q 81 (t).Because the solar flux, Q 81 (t), is a quasi-periodic function with a period of ∼ 11 years, for data sets longer than 11 years, the dominating competition is between the warming episode and trend.There is then a trade-off between the two best-fit values, which depend upon the observed values in the entire time series, i.e., data length.To see how this interdependence or correlation affects the ∼ 1 K decade −1 discrepancy more explicitly, we recall that the proxy of the warming episode is the function γ (z)P (z, t), which is shown in Fig. 3a for selected altitudes.This function rises to a peak temperature (max temperature), T p , at the time t p = t 0 + n(t 2 /t 1 )/[(1/t 1 )+(1/t 2 )].Comparing t p and T p in the 25-and 18-year-long data sets reveals the difference of their warming episode affecting the temperature trend.We plot these quantities as a function of altitude in Fig. 3b for F-11P(90-14_Avg) and for F-11P(90-07).Note that the altitude dependences for the two data sets are similar.Between 88 and 102 km, where the lidar signal is strong, we see little difference in t p but a systematic difference in T p between the two data sets.Above 93 km, t p is about constant, but T p increases continuously.Note that the peak warming (or maximum temperature response), T p , from the shorter data set is consistently higher by 0.5 to 2.5 K, depending on altitude, implying that a larger share of observed temperatures in the 1990s are attributed to the warming episode; this leads to a lower share for the trend assessment and thus a smaller cooling trend.With the longer data set, the reverse is true.Since the longer data set assesses the lingering warming episode more fully, it can render better judgment on the sharing between the two competitors; of course, more data leads to statistical accuracy and thus to a smaller uncertainty than the corresponding trend deduced from the shorter data set as shown in Fig. 2. We thus accept F-11P(90-14_Avg) from 25 years of Na lidar observations, marked in solid black circles, as the deduced midlatitude mesopause region temperature trend.

Discussions
The warming episode in our data plays a critical role in the temperature trend analysis based on our data sets.Furthermore, we assume a single temperature trend over 25 years in analysis.We thus shall discuss these issues before the final conclusion.

Mesopause warming and global surface cooling in the 1990s
A significant 6 K warming in 1992 and 1993 between 60 and 80 km was reported by Rayleigh lidar observations in southern France and attributed to the Mt Pinatubo eruption (Keckhut et al., 1995).Our suggestion (She et al., 1998) that the observed warming episode is one of the consequences of the Mt Pinatubo eruption does not yet have a clear geophysical causal relationship.To our knowledge, there has been no succinct explanation or model simulation published that re-lates the direct radiative and/or indirect dynamical effects of the Pinatubo eruption to the observed response in the mesosphere which lingered for ∼ 7 years at altitudes between 88 and 102 km as shown in Fig. 3a.There is, however, a comprehensive study of global mean surface temperature change in response to volcanic eruption and ENSO (El Niño-Southern Oscillation) events by Thompson et al. (2009).In response to the Mt Pinatubo eruption, they found a peak cooling of ∼ 0.3 K about 1.5 years after the Pinatubo eruption, which, remarkably, also lingered for ∼ 7 years.Changing the sign of their deduced global surface temperature response and multiplying it by a factor of 50, we obtain an episodic response very similar to our mesopause region temperatures response.
Better yet, we find that the scaled global surface temperature change, STR • (−50), deduced by Thompson et al. (2009) can be fitted to γ (z)P (z, t) that was used to represent the warming episode deduced from Na lidar observation.Next, we compare the peak delay time and the response time constant, t pd = t p −t MPE = t p −1.45, and τ = t 1 +t 2 , respectively, for the two observed episodes that occurred in the same time frame but ∼ 100 km apart in height.Thompson et al. (2009) analyzed the surface temperature response to a volcano T (t) by using the forcing function F (t) in terms of a simple model system of exponential decaying memory with time constant τ 0 = Cβ, where C is the effective heat capacity of the global atmosphericoceanic mixed layer per unit area and β is the damping coefficient, a measure of the climate sensitivity.They deduced C = 4.8 × 10 7 J m −2 K −1 , equivalent to the effect of the global atmosphere plus 9 m of the global oceanic mixed layer, and set β to be 2/3 K (W m −2 ) −1 leading to a time constant for the model system of τ 0 = Cβ = 3.2×10 7 s = 1.01 years.Though the system response to an impulse is a simple exponential decaying function with a memory of ∼ 1 year, the forcing function, with the Northern Hemisphere aerosol index as a proxy lasted several years through the end of 1994.The Aero Index (NH) • 100 is shown in Fig. 4a.Thompson et al. (2009) produced a cumulative response with a memory much longer than 1 year; the scaled response, STR • (−50), is also shown.Since STR • (−50) has a shape similar to the mesopause region warming episode, we fit it to the function γ (0)P (0, t) + background, giving γ = 13.0K, t 0 = 2.42 years, t 1 = 0.35 years, t 2 = 1.70 years, and t p = 2.88 years, along with a background of 0.39 K, shown in Fig. 4a as the blue curve, which is seen to match the scaled surface temperature response very well.The time constant and peak delay time for the global surface temperature response are, respectively, τ = 2.05 years and t pd = 1.43 years.Also in Fig. 4a is the warming response at 100 km in altitude (with the same background of 0.39 K added), the best-fit parameters of which are γ = 11.4K, t 0 = 2.81, t 1 = 0.20 years, t 2 = 1.60 years and t p = 3.17 years, or τ = 1.80 years and t pd = 1.65 years.It is evident that both blue and red curves have a similar shape, both spanning about 7 years, and thus can be represented by the same function.A more complete comparison between the warming episode in the mesopause region temperatures and the global surface temperature anomaly is shown in Fig. 4b, where the warming episode response time constant τ and peak delay time t pd are plotted as a function of altitude along with these values for the surface temperature response at the bottom of the figure.Averaged over the warming episode between 88 and 102 km, we find τ = 1.81 ± 0.42 years and t pd = 1.82 ± 0.26 years, compared with the surface temperature anomaly of τ = 2.05 years and t pd = 1.43 years.In addition, since the functional shape of both anomalies represents the distribution of the transit times of respective events, the concept of the "age of air" (AOA) (Waugh and Hall, 2002) is useful.Though the AOA concept was mostly applied to the transport of species from the tropical troposphere to the stratosphere, we use it here to describe the temporal history of an episodic response to a strong impulsive forcing.In fact, when area is normalized, the response of both surface temperature and mesopause temperatures after the Mt Pinatubo eruption is what is called the "age spectrum" in the AOA literature.Though all information on the transport process in question is contained in the age spectrum, the mean age (or the first moment of the spectrum in reference to the time of the forcing impulse, i.e., t MPE here) of air, t MA , may be used as a rough measure of the life of the process.With the "age spectrum" at each altitude, we can compute the associated mean age, t MA , also plotted in the figure, along with an upper-case M for the surface temperature response, t MA = 2.51 years.Averaged between 88 and 102 km, we have t MA = 2.92 ± 0.33 years for the warming response.
All these time constants are deduced from observational data.Assuming the Pinatubo aerosol reached the tropical lower stratosphere in negligible time as Mt Pinatubo erupted, for the warming episode, the mean age t MA is the time that the direct plus indirect (dynamic and feedback) effect of Pinatubo aerosol reaches the midlatitude mesopause region from the tropic lower stratosphere.For global surface cooling, the time t MA starts as the perturbation moves from the tropical lower stratosphere down through the global troposphere to the global atmospheric-oceanic mixed layer and the subsurface ocean (Thompson et al., 2009).We propose no detailed mechanism, which may be complicated, that leads to the mesopause warming after the Mt Pinatubo eruption.We hope that similarities in the observed and deduced response times between the volcanic eruption and the observed episodes, along with treating the surface temperature cooling with the same response function, will spur scientists and modelers to ascertain the causal effects of these strong episodic responses that occurred at the same time but ∼ 100 km apart in height.

A single linear trend or piecewise linear trends
The use of a single linear trend for a long data set is consistent with the classic recommendation of the World Meteorological Organization (WMO), using ∼ 30 years or more for analysis (Cnossen, 2012), and the practice of modelers, typically using 20 or ∼ 50 years of simulation for trend studies (Akmaev et al., 2006;Garcia et al., 2007;Berger and Lübken, 2011).It is nonetheless an assumption (Laštovička et al., 2012).Since a primary cause of a long-term temperature change is the anthropogenic emission of greenhouse gases, a single linear trend over a 25-year span may not reflect the reality of the rate of greenhouse emission changes.The emission of CO 2 and CH 4 into the atmosphere continues to increase.Indeed, Emmert et al. (2012) recently reported the observed increase of thermospheric CO 2 concentration.However, the loss rates (between 50 and 30 hPa) of the dramatic Antarctic ozone decrease (ozone hole, appeared in the late 1970s) have remained stable since 2000 (Hassler et al., 2011).In the midlatitudes, the trend began to reverse in the late 1990s (Akmaev, 2012;Qian et al., 2013), and it has been stabilized in recent years at a level below that in the 1960s.This recent O 3 rate change should slow down the cooling rate in the mesosphere and justify the use of a nonlinear (Keckhut et al., 2011) or piecewise linear trend for the regression analysis (Laštovička et al., 2012).Indeed, in analyzing the long-term variation in the reflection heights of radio waves from 1961 to 2009 (Bremer and Berger, 2002) and investigating temperature trends in the summer mesopause, Berger and Lübken (2011) and Lübken et al. (2013) found it appropriate to use three different linear trends with breaking points in 1979 and 1997.One could then reanalyze our 25year data using piecewise trends in the future.In this case, we would replace the term β(z)t in Eq. ( 1) by two different linear trends joined at a breaking point.Again, because the influences of solar flux, warming episode, and trends on temperatures are not independent, all these terms, along with the break point, if it is to be statistically determined, must be included in Eq. ( 1) to compete for the same nightly mean temperatures for regression analysis.

Conclusions
We have performed a regression analysis for the deduction of the mesopause region temperature trend based on an unprecedented Na lidar data set between March 1990 and December 2014.The 81-day averaged F10.7 solar is used as a proxy for solar activity, and a linear trend is assumed.Owing to a strong warming episode in the 1990s, the quarter century data set (90-14_Avg) is least-square fitted to an 11-parameter nonlinear model.The temperature trend shown in Fig. 2 starts from an insignificant value of 0.64 ± 0.99 K decade −1 at 85 km, increases to a maximum of 2.8 ± 0.58 K decade −1 between 91 and 93 km, and then gradually decreases to a warming trend above 103 km.Compared to the trend from the shorter data set, (90-07), which has a marginally significant cooling maximum of 1.55 ± 1.15 K decade −1 at 91 km (She et al., 2009), the quarter century data set gives a statistically quite significant cooling trend, larger by ∼ 1 K decade −1 .The discrepancy is due to the competition between warming episode and temperature trend.Since the longer data set assesses the longlasting warming episode more fully, it leads to more significant results.The mean uncertainty between 88 and 102 km are 0.6 and 1.3 K decade −1 , respectively, for the long and shorter data sets.The altitude dependence from the two data sets is quite similar, and their magnitudes are in the general range predicted by models.The trends reported here are −1.0 ± 1.0 K decade −1 at 87 km and −2.5 ± 0.5 K decade −1 at 90 km; they are consistent with the recently reported trends: respectively, −2.3 ± 0.6 K decade −1 (Offermann et al., 2010) and −4 ± 2 K decade −1 (Hall et al., 2012).
With regard to an interesting connection, we analyzed the surface temperature response after the Mt Pinatubo eruption reported by Thompson et al. (2009) with the same functional dependence as that used for the observed warming episode in the mesopause region.We determined the respective peak delay time, t pd , and mean age, t MA , to be, respectively, 1.43 and 2.51 years for surface temperature and 1.82 and 2.92 years for mesopause region warming.These similarities between the global surface temperature anomaly and the mesopause warming episode should hopefully spur community scientists and modelers to figure out the causal effects of these interesting phenomena.

Figure 1 .
Figure 1.Time series of nocturnal mean temperature recorded by a Na lidar at 86 km (a) and at 99 km (b).Included in (b) is also 81-day F10.7 solar flux in red with the times for Mount Pinatubo eruption (MPE), t MPE , and solar minima (Solar Min) marked.Data (black circles) between March 1990 and March 2010 were acquired at CSU (41 • N, 105 • W), and those (blue circles) between September 2010 and December 2014 were acquired at USU (42 • N, 112 • W).

Figure 2 .
Figure 2. Linear temperature trend from the quarter century data set with 11-and 7-parameter analyses, respectively denoted as F-11P(90-14_Avg) in black solid circles and F-7P(90-14_Avg) in black open circles.Shown for comparison are those data published based on an 18-year data set denoted as F-11P(90-07) in red solid squares and F-7P(90-07) in open red squares.

Figure 3 .
Figure 3. (a) Episodic warming in the 1990s at selected altitudes.By the beginning of 1999, the warming episode is basically over between 88 and 102 km.(b) The peak temperature (maximum temperature response), T p , and the time at which it occurs (or the time of max response), t p , of warming episode deduced from 25-year data set (90-14_Avg) is shown in black solid circles and open circles, respectively.Results shown in red solid squares and open squares were deduced from the 18-year data set (90-07).

Figure 4 .
Figure 4. (a) A comparison between 100 • stratospheric aerosol (black) and the scaled global surface temperature response, STR • (−50), in blue dots, from Thompson et al. (2009), along with its best fit (blue curve) and episodic warming response (EWR) in temperatures at 100 km (red curve).Both STR and EWR are over at the beginning of 1999, ∼ 7.5 years after the Mt Pinatubo eruption.(b) Deduced altitude-dependent time constant of the response, τ , peak delay time, t pd , and mean age, t MA , for the warming episodic function, respectively in open blue circles, open red circles, and solid black circles.Marked at the bottom of the figure are the times for the surface temperature response, respectively by blue and red crosses and the letter M. Data, STR • (−50) in (a) is derived from http://www.atmos.colostate.edu/~davet/ThompsonWallaceJonesKennedy/.