Lidar Observations of Instability and Estimates of Vertical Eddy Diffusivity Induced by Gravity Wave Breaking in the Arctic Mesosphere

On the night of 18–19 October 2018, sodium resonance lidar measurements show the presence of overturning in the mesospheric sodium layer. Two independent tracers, sodium mixing ratio and potential temperature, derived from resonance and Rayleigh lidar measurements, reveal that vertical spreading of the sodium mixing ratio contours and a layer of convective instability coincide with this overturning. Analysis of lidar measurements also reveals the presence of gravity waves that propagate upward, are saturated, and dissipate at the height of the convective instability. The vertical spreading is analyzed in terms of turbulent diffusive transport using a model based on material continuity of sodium. Estimates of the turbulent eddy diffusion coefficient, K, and energy dissipation rate, ε are derived from the transport model. The energy dissipated by the gravity waves is also calculated and found to be sufficient to generate the turbulence. We consider three other examples of overturning, instability and spreading on the nights of: 17–18 February 2009, 25–26 January 2015, and 8–9 October 2018. For all four events we find that the values of K (∼1,000 m2/s) are larger and the values of ε (∼10–100 mW/kg) are of similar magnitude to those values typically reported by ionization gauge measurements. These examples also reveal that higher levels of turbulent mixing are consistently found in regions of lower stability.


Introduction
Turbulence is important in understanding the coupling, energetics and dynamics through the whole atmosphere. In the mesosphere-lower-thermosphere (MLT) region, effects of turbulence include local heating as well as mixing and transport of heat, momentum, and constituents. Previous study has shown that turbulent transport is required to explain the enhanced vertical transport of nitric oxides from the thermosphere to the stratosphere in winter time (Smith et al., 2011). Whole atmosphere models have been used to assess the importance of diffusive transport and sensitivity analyses show that increasing the eddy diffusion coefficient by a factor of two yields significant changes in the vertical transport of minor species (Garcia et al., 2014;Meraner & Schmidt, 2016). The values of the eddy diffusion coefficient used in the whole atmosphere models vary between 10 and 100 m 2 /s in the MLT region, depending on altitude and latitude (Garcia et al., 2014;Meraner & Schmidt, 2016;Smith et al., 2011). Satellite observations of carbon dioxide and atomic oxygen in the MLT have yielded global-mean estimates of eddy diffusion coefficients that are consistent with these values (Salinas et al., 2016;Swenson et al., 2018). However, in situ rocket-borne measurements have shown a much larger range of turbulent eddy diffusion coefficients (between 1 and 1,000 m 2 /s) in the MLT, where the higher values occur in regions with instabilities (e.g., Bishop et al., 2004;Collins et al., 2011;Lehmacher et al., 2006Lehmacher et al., , 2011Lehmacher & Lübken, 1995;Lübken, 1997;Strelnikov et al., 2017;Szewczyk et al., 2013). Given that large variations exist in the measured values of turbulence under different environmental conditions, measurements of turbulent parameters under well-defined meteorological conditions are essential.
In the MLT region, turbulence is generated in instabilities that are generated by gravity waves (GWs) (e.g., Becker, 2012;Hines, 1988;Fritts & Alexander, 2003). As GWs propagate upward, wave amplitudes grow exponentially in the absence of dissipation as atmospheric density decreases with altitude. Eventually GWs become unstable and generate turbulence and deposit momentum and energy to the mean flow. However, the turbulence arising from GW dissipation remains one of the least quantified aspects of GW forcing in the MLT region (Fritts & Alexander, 2003). Collins et al. (2011, hereafter referred to as C11) used Rayleigh and resonance lidar measurements at Poker Flat Research Range, Chatanika, Alaska (65°N, 147°W) to study a strong local convective instability, which coincided with downward transport of sodium and increase of scale height of sodium. C11 analyzed the structure of the sodium layer in terms of vertical diffusive transport of sodium number density and estimated the turbulent eddy diffusion coefficient and energy dissipation rate, finding relatively large values of the eddy diffusion coefficient coinciding with convective instability. The approach of C11 differs from that of Gardner and Liu (2010) who report vertical transport due to fluxes associated with dissipating waves. These fluxes are directly calculated from measurements of sodium density and vertical wind perturbations and include the effects of all waves with vertical wavelength of 1-30 km and period of 3 min-15 h. Given the instrumental noise associated with the measurements, Gardner and Liu (2010) report average and seasonal profiles of the vertical fluxes based on averages of 370 h of observations acquired on 49 nights over 30 months.
In this study, we extend the scope of C11 by analyzing four sets of Rayleigh and resonance lidar observations on the nights of: 17-18 February 2009, 25-26 January 2015, 08-09 October 2018 October 2018 at PFRR. Our analysis includes the following four elements that did not appear in C11: First, we use potential temperature and sodium mixing ratio (rather than temperature and sodium density) to quantify the instability and spreading, respectively. Sodium mixing ratio is a more accurate tracer of material transport than sodium number density (Xu & Smith, 2003, 2005Xu et al., 2006). Potential temperature is an independent dynamic tracer; Second, we develop a turbulent transport model based on the mass continuity framework to describe the observed transport as turbulent diffusion and estimate the turbulent diffusion coefficient and energy dissipation rate; Third, we characterize GWs and quantify the saturation and dissipation of the waves. We compare the turbulence dissipation rates to the energy available from dissipating GWs and determine if energy from the GWs is sufficient to generate these levels of turbulence; Fourth, we draw on four case studies to investigate the relationship between the strength of turbulence and the atmospheric stability. In Section 2, we describe the Rayleigh lidar and resonance lidar systems and quantities measured with them. In Section 3, we present the observations of instability, vertical transport and GW activity on the night of 18-19 October 2018. In Section 4 we present our model of turbulent transport and estimate eddy diffusion coefficients and dissipation rates and compare them with the energy available from the GWs. In Section 5, we discuss these observations in terms of other studies and present our conclusions.

Rayleigh Density Temperature Lidar
The Rayleigh Density Temperature Lidar (RDTL) was installed at PFRR in 1997. The transmitter of the RDTL system incorporates a Nd:YAG laser operating at 532 nm and has remained the same. The receiver has been upgraded from a single channel system (Collins et al., 2011), to a two-channel system in 2015 (Triplett et al., 2018), and a three-channel system in 2017 (Alspach, 2020). The raw lidar profiles were acquired at a temporal and spatial resolution of 50 s and 75 m for the single-channel receiver and 50 s and 48 m for the dual and triple channel receivers. The raw data are integrated in time and smoothed in altitude to increase the signal-to-noise ratio in the density and temperature retrievals. Different intervals are selected to retrieve different parameters. We retrieve temperature profiles over 2 h intervals at 15 min steps and use these profiles to characterize the temperature structure in the stratosphere and mesosphere and calculate potential temperature profiles. Under the assumption of hydrostatic equilibrium, we derive the potential temperature from the temperature profile (Franke & Collins, 2003). We calculate atmospheric density profiles over 1 h intervals at 15 min steps, which is sufficient to obtain robust relative density profiles up to 90 km. The relative density profiles are normalized by radiosonde density measurements over the 30-32 km altitude range to obtain absolute density profiles. To characterize gravity-wave activity in the stratosphere and mesosphere, we choose a higher temporal resolution and calculate density profiles over 30 min intervals at 5 min steps in order to include higher frequencies of the GW spectrum.

Resonance Lidar
We use measurements from different sodium resonance lidar systems in this study. The measurements in 2009 and 2015 were made with a broadband sodium density resonance lidar (Triplett et al., 2018), while the measurements in 2018 were made with a narrowband sodium resonance wind-temperature lidar .

Broadband Sodium Resonance Density Lidar (SRDL)
SRDL observations were made with an excimer-pumped dye laser since 1999. The lidar measurements yield sodium number density profiles (∼70-120 km) using established techniques (Collins et al., 2011). The raw lidar profiles were acquired at a resolution of 50 s and 75 m in 2009 and 100 s and 75 m in 2015. In this study we calculate sodium profiles over 1 h interval at 15 min steps. We then take the ratio of these sodium number density profiles to the absolute atmospheric density measured by the RDTL to calculate the sodium mixing ratio. The sodium mixing ratio profiles are therefore obtained in the altitude range of 70-90 km, where both Rayleigh and resonance lidar data overlap in altitude.

Sodium Resonance Wind-Temperature Lidar
The Sodium Resonance Wind-Temperature Lidar (SRWTL) was installed at PFRR in 2017 and ongoing measurements of the sodium number density, temperature and meridional wind in the MLT (∼70-120 km) have been acquired since then . A zonal wind channel was added in 2019 after the observations presented in this paper. The SRWTL is a narrowband sodium resonance lidar that employs a three-frequency technique and a pulse-dye amplifier architecture to measure the wind and temperature (Bills et al., 1991;She et al., 1990She et al., , 2002She & Yu, 1994. The light from the transmitter is split into two beams with one pointing zenith and the other pointing 20° off-zenith toward the north. The receiver consists of two telescopes pointing in the same directions as the transmitted beams. This system provides measurements of sodium number density, neutral temperature, meridional wind and vertical wind. The raw lidar profiles were acquired at 8 s and 30 m. We integrate the signal over 5 min intervals and spatially smooth the data over 500 m in the vertical to reduce the statistical uncertainty and determine temperature and wind. The retrieval yields temperature and wind profiles with a resolution of 5 min and 500 m which are used to characterize GWs. In addition, the temperature data integrated over 2 h intervals at 15 min steps to provide the initial temperature estimate at 87.5 km for the RDTL temperature profiles (Section 2.2.1). The high resolution and accuracy of the SRWTL temperature measurements allow us to characterize atmospheric instability at higher resolution than previous studies where Rayleigh lidar temperature measurements were used to characterize instability and wave activity (e.g., Triplett et al., 2018).
The retrieval of the sodium mixing ratio uses measurements from both SRWTL and RDTL systems. The SR-WTL temperature measurements in the MLT region (∼70-120 km) allow us the extend the RDTL absolute density retrieval, and thus the sodium mixing retrieval, into this region. To achieve this extension, we first combine the temperature profiles measured by these two lidars to form continuous temperature profiles from 30 to 120 km at 1 h intervals and 15 min steps. These temperature profiles are then used to calculate the density of the atmosphere under the assumption of hydrostatic equilibrium. The lidar measured density is then normalized by the density measured by radiosonde to yield the atmospheric density profile. The sodium mixing ratio (15 min step, 1 h interval) is retrieved as the ratio of sodium number density to absolute atmospheric density in the altitude range of 70-105 km.

Sodium, Temperature, and Wind
We use the observations on the night of 18-19 October 2018 to illustrate the overturning, instability, spreading, and wave dissipation. The sodium number density, temperature, and meridional wind measured by the SRWTL on this night are shown in Figure 1. The sodium number density and temperature are measured by the vertical channel of the SRWTL, while the meridional wind is derived from the combination of the vertical and meridional channels. There is an overturning signature (i.e., lower density on top of higher density) on the bottom side of the sodium layer around 20:15 LST and 87 km (black box in Figure 1a). Meanwhile, there is a large negative gradient (∼-13 K/km) in the temperature ( Figure 1b) and a strong shear (-40 m/s/ km) in the meridional wind ( Figure 1c) which coincide with the overturning in the sodium number density and indicate the presence of reduced convective and dynamic stability in this region. Indeed, Figure 1d shows the broadening and steepening of the potential temperature contours at 20:15 LST near 87 km (the black square box). The sodium mixing ratio plotted in false color (Figure 1d) follows the potential temperature and shows broadening and steepening that confirms the behavior of the potential temperatures at 20:15 LST and 87 km. The steepening of the potential temperature contours (i.e., contour lines become vertical) indicates reduced convective stability (Franke & Collins, 2003), while the broadening of the sodium mixing ratio contours (for instance, the distance between the 5 parts per trillion (PPT) and 20 PPT contour lines increase from a value of 2 km before 20:00-6.5 km at 20:45 LST) indicates vertical transport and mixing. The fact that the potential temperature and sodium mixing ratio show consistent behavior confirms that the convective instability coincides with material transport.
Signatures of upward propagating GWs with relatively continuous downward phase progression are found below ∼87 km and disappear above 90 km around 20:15 LST (Figures 1b and 1c, black box), which implies that the waves are breaking or being strongly dissipated at these altitudes and time. These wave-breaking signatures are more carefully examined in Section 3.3, where the perturbations due to the wave are investigated. The coincidence of the overturning signature, convective instability, transport and changes in GW (highlighted by the black box in Figure 1), suggests that all these phenomena are related to each other.

Instability and Overturning
We show the vertical profile of the sodium mixing ratio (solid, yellow) and sodium number density (dotted, yellow) averaged over 60 min centered at 20:15 LST in Figure 2a. There is a well-defined 3.5 km deep region between 86.5 km and 90 km where the sodium mixing ratio spreads out (i.e., the reduction in vertical gradient in solid orange line) and the sodium number density overturns (i.e., lower density above higher density of the orange dotted line around 87 km). We see that the region of spreading and overturning corresponds to the region of reduced convective stability in Figure 2b, which shows the vertical gradient of the potential temperature (  / d dz), calculated over 5,15, 30, and 60 min intervals centered around 20:15 LST. The region of sodium density overturning and sodium mixing ratio diffusion is broader than and spans the region of LI ET AL. convective instability. The gradient is calculated as the slope of a linear fit to the potential temperature profile over a 2 km interval at 30 m steps. The reduction of stability around 87 km is evident at all temporal scales which implies that the reduction of stability in this layer is a robust feature. It is the most pronounced in the highest resolution 5 min profile (i.e., most unstable). The instability region in the 5 min profile where  / d dz < 0 is 0.8 km deep extending from 87.2 to 88.0 km. Modeling studies have shown that this behavior in the sodium mixing ratio and number density is consistent with waves that are approaching or have reached instability (Xu et al., 2006).

Gravity Wave Activity
We illustrate the gravity wave activity with the temperature perturbations from the SRWTL measurements (80-105 km) and the relative density perturbations from the RDTL measurements (65-90 km) in Figure 3. The temperature and density perturbations are calculated as the difference between the temperature and density profiles and their nightly average profiles respectively, and the relative density perturbations are calculated as the ratio of the density perturbations to the nightly average density profile. The temperature and relative density perturbations show similar wave periods with clear downward phase progressions corresponding to upward-propagating waves from ∼75 to ∼90 km. The temperature perturbations in the MLT range from +25 K to −25 K ( Figure 3a) and the atmospheric density perturbations in the upper stratosphere and mesosphere reach amplitudes of 6% ( Figure 3b). Consistent with GW polarization relation, the temperature and density perturbations appear out of phase, with the warm (cold) phase of the temperature perturbations associated with the less (more) dense phase of the density perturbations.
To identify the dominant GWs, we determine the frequency spectra for sodium number density, temperature and meridional wind perturbations at each altitude. The frequency spectra power is normalized by the total spectral power at each altitude to produce the normalized spectra power (Figure 4). There is a 5.6 h wave between 80 and 90 km in all three spectra. This 5.6 h wave disappears around 90 km, which coincides with the location of the overturning/instability structure (Figure 1), implying that the wave breaks around this altitude. An interesting feature is that in the region where the 5.6 h wave disappears, a 10 h wave appears in all three spectra between 90 and 94 km. Given that the period of this wave is approximately twice of the 5.6 h, this feature may imply that the 5.6 h wave only breaks in certain phases of the wave, resulting LI ET AL. in a doubling the observed period. For instance, in Figure 3a, between 90 and 93 km, the warm temperature phases in the beginning and at the end of the night tend to follow the phase progressions from lower altitudes, while the warm phase in the middle of the night is smeared out and replaced by cold phase. The phenomena of GW breaking appearing in one phase of the wave is consistent with recent model results (Pestana et al., 2018).
We determine the amplitude and phase of this wave by fitting a 5.6 h harmonic to the perturbations in sodium relative density, temperature, meridional wind, and atmospheric relative density at each altitude (Figure 5). The sodium relative density is derived by dividing the sodium density by the nightly mean sodium density profile. The wave amplitude grows exponentially with altitude in the upper mesosphere, maximizes at ∼87 km, and then decreases with altitude. The wave phase shows a downward phase progression below ∼90 km, and upward phase progression above. Such transitions in wave amplitude and phase suggest that the ∼5.6 h GW experiences significant dissipation and likely breaks near ∼90 km. The vertical wavelength derived from the vertical phase progression is 11.6 km from RDTL density, 11.6 km from SRWTL sodium relative density, 10.4 km from SRWTL temperature, and 11.5 km from SRWTL meridional wind in the range of 80-90 km, where the downward phase progression corresponding to upward energy propagation is observed. We summarize the observed GW characteristics in Table 1. The observed period and vertical wavelength indicate that this is an inertia-GW rather than a tidal wave because the vertical wavelengths of subdiurnal tides are much longer. Theoretical calculations and model simulations by Smith and Ortland (2011) indicate that 8 h tides have vertical wavelength larger than 47 km, while observations and simulations by Smith et al. (2004) reveal that shorter period tides (6 h, 4.8 h, etc.) have even longer vertical wavelengths in the Arctic MLT region. Furthermore, the amplitudes of these tides are much smaller than the measured wave amplitudes (∼1 K in temperature, 2 m/s in meridional wind).
We determine the wave amplitude in horizontal wind using the GW polarization and dispersion relations. The equations are introduced in the Supporting Information S1. We assume that the ground-based wave frequency, ω, is nearly equal to the intrinsic frequency,  . The meridional wind (  v) amplitude determined from the SRWTL measurements is 19.2 m/s. From Equation S1.8, we calculate the wave amplitude in the horizontal wind is  H u is 32.6 m/s and thus the zonal wind (  u) amplitude is 26.3 m/s. The horizontal wavelength, λ h , of this wave is calculated as 777 km (Equation S1 .3) and the horizontal phase speed, c H , is calculated as 39.5 m/s (Equation S1.9). The GW characteristics indicate that this wave is likely approaching linear instability as the ratio of the wave amplitude ( 32.6 m/s H u   ) to the phase velocity (c H = 39.5 m/s) is over 80% (e.g., Staquet & Scmmeria, 2002). This suggests that the 5.6 h wave propagates upward to the mesosphere and breaks near 90 km, which is consistent with the disappearance of the wave in the spectrum (Figure 4) and the vertical structure of wave amplitude and phase ( Figure 5). Since  H u derived from the ground-based frequency (ω) points to the existence of wave saturation and instability, using ω to approximate  appears to be a reasonable assumption.
We calculate the potential energy per unit mass or specific potential energy (PE) for both the monochromatic wave and the ensemble of waves using the temperature perturbations measured by the SRWTL at 5 min resolution (e.g., Lu et al., 2009;Thurairajah, Collins, Harvey, Lieberman, Gerding, et al., 2010;Triplett et al., 2018). Using the amplitude of the 5.6 h harmonic fit to represent the temperature perturbation term, the PE of the monochromatic GW decreases to near zero (∼1-2 J/kg) around 93 km ( Figure 6). For waves of these periods, the kinetic energy (KE) of the wave is nearly equal to the PE (e.g., Nappo, 2002;Geller & Gong, 2010). Given the average PE of the wave over 86.5-90.0 km is 260 J/kg, the total energy (i.e., PE + KE) of the wave is 520 J/kg. To derive the PE of the ensemble of waves, we first apply a two-dimensional Fast Fourier Transform based filter in time and altitude to the temperature perturbations (Figure 3a) to remove spectral components with periods longer than 11 h and vertical wavelengths longer than 15 km from the ensemble of waves. The filtering also removes the downward propagating GWs and retains the upward propagating waves (e.g., Lu et al., 2017Lu et al., , 2015aZhao et al., 2017). We calculate the PE for the ensemble of waves from the filtered temperature perturbations (e.g., Chu et al., 2018;. Like the monochromatic wave, the PE of the ensemble of waves decreases with altitude in the altitude range from 88 to 92 km (blue line in Figure 6). The altitude ranges where the PEs of both the monochromatic and ensemble of waves decrease overlap the region of overturning in the sodium layer from 86.5 to 90 km (Figure 2). The PE at the lower altitude of the layer of sodium spreading (i.e., 86.5 km) is 330 J/kg and the PE at the upper altitude (i.e., 90 km) is 274 J/ kg. Thus the total energy of the waves is 660 and 548 J/kg at the lower and upper altitude respectively. These numbers allow us to quantify the energy dissipated by the GWs in Section 4.2.

Estimation of Eddy Diffusion Coefficient
Turbulence has been observed coincident with the presence of convectively unstable layer (e.g., Lehmacher et al., 2006;Lehmacher & Lübken, 1995;Lehmacher et al., 2011;Szewczyk et al., 2013;Triplett et al., 2018). We develop a turbulence transport model based on the continuity equation of sodium number density, which is expressed as (Jacob, 1999), where n s is the sodium number density and the three terms on the righthand side (RHS) represent the time changing rates of sodium number density induced by advection (adv), turbulence (turb), and chemistry (chem), respectively. These terms can be derived as where n a is atmospheric number density, f = n s /n a is sodium mixing ratio,  K is matrix form (tensor) of eddy diffusion coefficients, P is chemical production, and L is chemical loss. We consider the one-dimensional situation based on transport in the vertical direction (z) alone. We assume that there is no horizontal transport will be discussed in Section 4. In this case, the turbulence term can be expressed in terms of the vertical gradient of the vertical turbulent flux (ϕ), LI ET AL.
where H f is the scale height of the sodium mixing ratio, and K is the vertical turbulent diffusion coefficient. The vertical turbulent flux of the sodium number density or eddy sodium flux is given by, We assume that the atmosphere was at a steady state before turbulence is induced by the wave breaking. The steady-state density of the atmosphere, sodium atoms, and sodium mixing ratio vary exponentially with height, where n s0 and n a0 are the number density of sodium and atmosphere at z 0 , respectively, and First, we assume that the advection term is negligible (the uncertainty from this assumption is assessed in Section 4.1), then Equation 1 becomes: s s a n f n Kn t z z t Assume that the eddy diffusion coefficient is K 0 and sodium number density is s n before instability-driven turbulence is induced. Then Equation 10 or the steady state can be written as: s s a n f n K n t z z t Now consider the case where extra turbulence is initiated at a certain time, t 0 , and an additional eddy diffusion coefficient, K', is induced. The perturbation induced by the turbulence is given by, n s '. Thus    The chemistry term can be expressed as where τ is the chemical time constant. Notice that n s ' is initially induced by turbulence thus it has the same , which can be positive or negative, depending on whether sodium is transported from a region of higher to lower or lower to higher mixing ratio by turbulence. The minus sign in Equation 14 indicates that the sign of is opposite to that of n s ', and indicates that chemistry counteracts the perturbation induced by turbulence. Initially the amplitude of the perturbation, n s ', is small, thus the chemical term, which is directly proportional to n s ' (Equation 11), is smaller than the turbulent term, which is dominated by the vertical gradient of s n (first term on the RHS of Equation 13). As time goes, the amplitude of n s ' increases with time, which causes the amplitude of the chemical term to increase. The effect of turbulence is to diffuse the sodium atoms, thus reducing the gradient of the sodium mixing ratio, thereby from Equation 5 the amplitude of the turbulent term decreases with time. Eventually, at a certain time (t m ), the turbulent and the chemical terms in Equation 13 balance, and n s ' reaches a local maximum.
To summarize, our model of turbulent diffusive transport of sodium simulates the dynamical process that turbulence triggers the change of n s ' initially and the amplitude of n s ' increases until it reaches a maximum at t m . During this period, turbulence causes vertical diffusion of the sodium atoms and thus increases the scale height of sodium mixing ratios. The chemical effect which is small initially increases with time as |n s '| increases and ultimately balances the turbulence effect at t m . After t m , |n s '| decreases because the turbulence effect is now smaller than the chemical effect, while the scale height keeps increasing until the turbulence effect disappears. Thus there is a delay between the appearance of local maximum of |n s ' | and the maximum of H f . Notice that n s ' can be both positive or negative, thus the value of n s reaches a local extremum at t m . We use these signature (i.e., the simultaneous increase of |n s ' | and H f , and the interval between their maxima) to identify the turbulent event and to find the time, t m , when the when the turbulence and chemical effects balance. This balance point allows us to estimate the eddy diffusion coefficient as follows. First,    / s n t becomes zero in Equation 13 at t m , combining with Equation 14, we have, Following previous studies, we assume a single slab of uniform turbulence over a layer of depth L, where the turbulent flux is confined within the layer and decreases to zero at the bottom of the layer (e.g., Collins et al., 2011;Whiteway et al., 1995). We integrate Equation 15 over the turbulent layer (of depth is L) and assume that at t m ,    s s n n. Then the expression for the eddy diffusion coefficient becomes, where the perturbation ratio, α, is defined as, Under the assumption that the turbulent term is negligible after t m , the decreasing rate of sodium mixing ratio after t m gives the chemical time constant τ.
LI ET AL.

10.1029/2020JD033450
We can now use Equation 16 to calculate the eddy diffusion coefficient. First, we determine the time of maximum (t m ), from the variation of the sodium mixing ratio with time in the body of the spreading event (Figure 7, orange solid line). We calculate the sodium mixing ratio averaged between 85 and 90 km (Figure 7, orange line) to capture the whole overturning structure (Figure 2, highlighted region, 86.5-90 km), and to minimize the GW influence on the calculation of n s ' . The scale height of the mixing ratio is calculated over the range of 86.5-90 km (Figure 7, black dashed line). The significant increase in scale height is one of the signatures of the event. Hence, we pick the period between 19 and 23 UT as the interval when the turbulence event occurs. The average sodium mixing ratio reaches a local maximum of 30.9 PPT at 20:15 LST, which we choose as t m . The scale height of sodium mixing ratio has a typical value of 3 km before the spreading event begins (19:30 LST), and then significantly increases to 8.1 km at 20:15 LST and reaches a local maximum of 11.6 km at 20:45 LST.
We assume that the chemical term dominates the turbulent term during the period immediately following t m from 20:15 LST to 22:45 LST. We calculate the chemical time constant, τ, from the decrease in the sodium mixing ratio by fitting a decaying exponential to the mixing ratio over this period. Notice that this mixing ratio is averaged over a 5 km vertical window as aforementioned. We find a value of τ of 75 min. We have also calculated the value of τ from the mixing ratio averaged over a 3.5 km range and the value is 76 min. The fact that the two values are so close and the fact that the vertical wavelength of the wave is ∼10 km suggest that by averaging over 5 km, the transport induced by the GW is significantly suppressed. At t m , the scale height, of sodium number density H s , atmospheric density, H a , and the sodium mixing ratio, H f , are −94.0, 6.9, and 8.1 km respectively. The corresponding mean value of n s is 3,847 cm −3 . The value of s n , calculated using the average of n s over the whole observation period, is 3,135 cm −3 . The corresponding value of α is 0.19. The depth of the turbulent layer, L, is also required to calculate the value of K'. As we saw in Section 3.1.1 (Figure 2), the depth of the overturning layer in sodium density is 3.5 km, while the depth of the convective instability is 0.8 km. We tabulate these values in Table 2. Using Equation 16, we estimate an eddy diffusion coefficient of 2.6 × 10 2 and 1.2 × 10 3 m 2 s −1 , corresponding to values of L of 0.78 and 3.5 km, respectively. These values of K' are tabulated in Table 3. The corresponding values of the eddy sodium flux are 1.0 × 10 8 and 4.6 × 10 8 atom m −2 s −1 .
Following Weinstock (1981) we calculate the turbulent energy dissipation rate (ε) from the eddy diffusion coefficient, K', as, where N is the Buoyancy frequency calculated as, where θ is the potential temperature, g is the gravitational constant, T is the temperature, Γ is the dry adiabatic lapse rate (9.5 K/km), and γ is the lapse rate of temperature profile (e.g., Dutton, 1986). We use temperature profiles with resolution of 5 min to estimate γ as 9.3 K/km (Table 2). This estimate of γ is calculated over a 3 km range from 86.1 to 89.1 km centered at the unstable layer (Figure 2, blue line) to avoid a negative N 2 due to a super adiabatic lapse rate at the center of the convective instability. N 2 is calculated to have a value of 1.0 × 10 −5 s −2 , which yields an estimate of ε of 2.1 (9.4) mW/kg using both estimates of K' (Table 3).

Gravity Waves as Source of Energy
To determine the energy source for the turbulence, we investigate the total energy associated with both the monochromatic GWs and the ensemble waves. In the following analysis, we assume that the diminishment   of the wave energy with altitude is totally attributed to wave dissipation. Certainly, there are other possible mechanisms that may cause the loss of the wave energy with altitude (e.g., evanescence and/or reflection of waves, waves propagating out of the measurement volume). However, this assumption is sufficient for investigating if the energy from the waves is enough to generate the observed level of turbulence. From the GW observations presented in Section 3.3 (Figures 5 and 6) we see that both the monochromatic wave and the ensemble of waves are not propagating freely with altitude but are losing energy as they propagate. We estimate the wave energy available for generation of turbulence by assuming that the specific potential energy of the wave is dissipated in the time it takes for the wave to travel through the layer. This dissipation time, τ d , is given by the ratio of the depth of the layer, L d , to the vertical group velocity, c gz . If the GWs loose an amount of energy, ΔE, then the GW energy dissipation rate is given by, The monochromatic wave is nearly completely dissipated at 93 km where we assume that the wave breaks ( Figure 6). Using the monochromatic GW parameters (total energy of 521.2 J/kg, c gz of 0.46 m/s, Table 4) and the depth of the dissipation layer layer, L d , of 3.5 km (Table 2), we calculate that the dissipation time is 127 min and the energy dissipation rate associated with wave breaking is 70 mW/kg. This value is larger than the value of 2.1 (9.4) mW/kg that we have estimated for the turbulent energy dissipation rate. The ensemble of waves also shows a decrease in PE from 88 to 92 km ( Figure 6). The total energy at the lower edge of the layer is 666 J/kg and at the upper edge of the layer is 548 J/kg. If the waves propagated freely, at the upper edge they should have an energy of, where E 0 is the PE of GWs at the lower edge, z 0 . Then the amount of energy dissipated, ΔE, is given by where E obs is the PE of GWs at the upper edge, z. For an atmospheric scale height, H a , of 7 km, we calculate ΔE = 545 J/kg. The representative ensemble of GWs has a spectrally weighted average period and vertical wavelength of 4.5 h and 8.7 km, respectively, and a corresponding vertical group velocity of 0.47 m/s (Equation S1.10 in Supporting Information). Given a dissipation depth, L d , of 3.5 km, we estimate a dissipation time, τ d , of 124 min and an energy dissipation rate, ε GW of 74 W/kg. This value of ε GW is again larger than the value of 2.1 (9.4) mW/kg that we have estimated for the turbulence dissipation energy rate. Therefore, the GW dissipation rate is sufficient to generate this level of turbulence. We summarize the characteristics of the monochromatic and ensemble GWs in Table 4. The fact that the GW energy dissipation rates are larger than the turbulent energy dissipation rates implies that dissipated GWs may generate a variety of other motions in addition to turbulence (e.g., secondary waves, Bossert et al., 2017;Becker & Vadas, 2018).

Three Other Case Studies
To investigate the relationship between strength of turbulence and environmental instability, we analyze the other three cases that also show the signatures of overturning and reduced stability on 08-09 October 2018, 25-26 January 2015, and 17-18 February 2009 (circled by the black boxes in Figure 8). The potential temperature for the 08-09 October 2018 case is derived from the temperature measured by the SRWTL, while the potential temperature for the 2019 and 2015 cases are from the temperatures measured by the RDTL. In Figure 9, we show the vertical profiles of dθ/dz, sodium mixing ratio and sodium number density for each case. Similar to the 18-19 October 2018 case, there is overturning in the sodium density, spreading in the mixing ratio and well-defined layers of convective instability. The values of the depth of the unstable layers varies between 0.5 and 2 km, while the depth of the layers of overturning varies between 2 and 3.5 km. Again, the depth of the overturning and spreading is larger than the depth of the instability. The values of the sodium time constants are between 52 and 79 min consistent with the chemical time constants for the sodium layer (Collins et al., 2011;Xu & Smith, 2003, 2005. The values of the perturbation ratio, α, are between 0.13 and 0.55. The lower values of α are found closer to the peak of the sodium layer (∼89 km) where the undisturbed mixing ratios are higher, and the highest value is found on 17-18 February 2009 when sodium is transported downwards near the lower edge of the sodium layer (∼77 km). We tabulate these values in Table 2.
We find values of the eddy diffusion coefficients between 140 and 920 m 2 /s based on the depth of the convective instability and between 390 and 1,600 m 2 /s based on the depth of the overturning in sodium density.
The corresponding values of the eddy sodium flux vary between 3.3 × 10 7 and 1.1 × 10 8 atom m −2 s −1 based on the depth of the instability and between 1.1 × 10 8 and 4.6 × 10 8 atom m −2 s −1 based on the depth of the overturning. The corresponding energy dissipation rates vary between 2.1 and 38 mW/kg based on the depth of the instability and between 7.0 and 108 mW/kg based on the depth of the overturning. The largest values of the eddy diffusion coefficient are found on 25-26 January 2015, when the sodium scale height is largest, the depth of the instability and the overturning is greatest, and the square of the Buoyancy frequency is smallest.
On 8-9 October 2018 we also characterize the GW activity in the MLT with measurements from the SRWTL (Table 4). We find the energy dissipated by the monochromatic GW is 83 mW/kg while the energy from the ensemble of GWs is 6.6 mW/kg. The energy deposited from the LI ET AL. The values outside the brackets are calculated using the depth of the unstable layers, while those inside the brackets are using the depth of the layers of diffusion.

Comparison with Other Observations
We plot the estimated turbulent eddy diffusion coefficients and energy dissipation rates along with multiple rocket-borne ion-gauge turbulent investigations in the Arctic in Figure 10. We consider four studies (Lübken, 1997  The L97-W profile represents the average of 12 wintertime measurements over two winters at ASC. The L97-S profile represents the average of 7 summertime measurements over three summers at ASC. The values obtained using the depth of the instability layers are plotted as red circles (Letal20 1), while the values obtained using the depth of the overturning of sodium density are plotted as blue circles (Letal20 2). For the estimates of turbulence using the depths of the instabilities, the values of K' (Figure 10a, red circles) are larger than the rocket-measured values while the values of ε (Figure 10b, red circles) are in good agreement with the values reported by ionization gauge measurements at PFRR. These comparisons suggest that it is more accurate to use the depth of the instability (0.5-2.0 km) as the depth of the turbulent layer (L in Equation 16) than the depth of the overturning of sodium density. This is consistent with the model results (e.g., Fritts et al., 2016.) and rocket-borne measurements (e.g., Triplett et al., 2018), which reveal that the depths of individual turbulent layers are several hundred meters which are similar to the depth of the unstable layers ( Table 2). The values of the eddy sodium flux are similar in magnitude to the values of the fluxes derived directly from wave-driven fluctuations of sodium density and vertical wind (10 7 -10 8 atom m −2 s −2 , Gardner & Liu, 2010). However, we note that the eddy values reported here are associated with individual case studies of wave breaking, while those reported by Gardner and Liu (2010) represent average values over 49 nights of observation.
For the four cases, the values of K' vary substantially. The calculation of K' is independent of N 2 , and for the three cases that are measured at similar height (∼85-90 km), larger K' corresponds to smaller N 2 (Table 3). This implies that a less stable atmosphere is more likely to generate stronger turbulent diffusion. The case on 17-18 February 2009 at a lower height (∼75 km) has much stronger turbulence with both K' and ε having much larger values than the other studies. While L97-W represents the average of 12 wintertime profiles, one of the 12 profiles shows a level of turbulence at 75 km that is similar to the values determined in this study ( Figure 6 of Lübken, 1997).
The estimates of the turbulent diffusion coefficients assume that the observed loss time of the sodium mixing ratio represents the chemical time constant, τ. This estimate is made under the assumption that advection is negligible. The observed lifetime may be significantly different from the true lifetime of the sodium if the sodium is advected by the background wind. For a layer with horizontal scale, L h , and vertical scale H v , in the presence of horizontal velocity U h , and vertical wind, U z , the observed lifetime, τ obs , is related to the true lifetime, τ, as   wavenumber, λ h and λ z . We use the case on 18-19 October 2018 as an example to quantify the uncertainty on the estimate of τ induced by ignoring the advection term. In the region of wave-breaking, the SRWTL measures meridional wind of 40 m/s. If we assume a zonal wind of the same magnitude, the horizontal wind speed would be 57 m/s. Thus we use typical horizontal wind speed of 50 m/s, which is consistent with other observations (e.g., Larsen et al., 2003;Lehmacher et al., 2011;Mesquita et al., 2020) and horizontal wavelength of 777 km (see Section 3.3), we estimate a value of τ h of 41 min. Similarly, using the corresponding wave vertical wind amplitude of 0.5 m/s and vertical wavelength of 10 km, we estimate a value of τ v of 53 min. With Equation 23, we calculate the true lifetime between 18 and 126 min. Thus we note that there is an uncertainty of a factor of four in our estimate of, τ v , and thus in our estimate of the turbulence diffusion coefficient. However, the values of τ from our measurements are consistent with the time constants determined from chemical models of the sodium layer (Xu & Smith, 2003, 2005, and thus the conclude that the actual uncertainty is less.

Summary and Conclusions
We have investigated overturning, convective instability, wave dissipation, and turbulence in the MLT based on lidar measurements of sodium, temperature, and density. We present results from four case studies based on nights of observations at PFRR, Chatanika, Alaska. We find that overturning in the sodium density coincides with instability and vertical spreading of sodium and potential temperature. The instabilities are embedded in the overturning layers in sodium density where the depth of the unstable layers (0.7-1.5 km) are less than the depth of overturning layers (2.0-5.0 km) in sodium density. We find that the overturning events coincide with dissipating GWs, and the energy dissipated by monochromatic waves is sufficient to generate the turbulence. We have developed a turbulent transport model based on material continuity of the sodium mixing ratio to explain the vertical diffusion. Based on this model we have estimated values of the turbulent eddy diffusion coefficient, K' and energy dissipation rates ε. The values of these coefficients are in better agreement with values derived from rocket-borne measurements when we use the smaller depth of the unstable layers instead of the depth of overturning layers in sodium density as the depth of the turbulent layers. The values of K' (140-920 m 2 /s) are larger than typically reported while the values of ε (2-38 mW/kg) are in good agreement with the values reported by ionization gauge measurements at PFRR. We find that larger values of K' correspond to lower values of N 2 and greater spreading of sodium. These relatively large values of the diffusion coefficients reflect the enhanced mixing due to convective instability and wave breaking. Figure 10. Summary of turbulent eddy diffusion coefficients (a) and energy dissipation rate (b) measured in the Arctic MLT region (see text for details). Values derived in this study are marked by red circles (derived using the depths of instability as the depths of turbulent layers) and blue circles (derived using the overturning depths as the depths of turbulent layers). Legends L97-W, L97-S, Letal11, Setal 13, Tetal18, Tetal18 Avg, and Letal20 O, denote for results from winter-time measurements from Lubken et al (1997), summer time measurements from Lübken et al (1997), Lehmacher et al. (2011), Szewczyk et al. (2013, Triplett et al. (2018), averaged measurements from Triplett et al (2018), the measurements in this study using the depth of the layer of overturning as the depth of turbulent layer; Letal 20 I, the same as Letal20 O except using the depth of the unstable layer as the depth of turbulent layer. MLT, mesosphere-lower-thermosphere.