Retrieval of atmospheric static stability from MST radar return signal power

An empirical technique for retrieving profiles of the square of the Brunt-V̈ ais̈alä frequency,ω2 B , from MST radar return signal power is presented. The validity of the technique, which is applied over the altitude range 1.0–15.7 km, is limited to those altitudes at which the humidity contributions to the mean vertical gradient of generalised potential refractive index, M, can be ignored. Although this is commonly assumed to be the case above the first few kilometres of the atmosphere, it is shown that humidity contributions can be significant right up to the tropopause level. In specific circumstances, however, the technique is valid over large sections of the troposphere. Comparisons of radarand (balloon-borne) radiosonde-derived ω2 B profiles are typically quantitatively and qualitatively well matched. However, the horizontal separation between the radar and the radiosondes (which were launched at the radar site) increases with increasing altitude. Under conditions of mountain wave activity, which can be highly localised, large discrepancies can occur at lower-stratospheric altitudes. This demonstrates the fact that radiosonde observations cannot necessarily be assumed to be representative of the atmosphere above the launch site.


Introduction
The vertical gradient of potential temperature gives a measure of the atmosphere's static stability, i.e. of its resistance to vertical motions.It is conveniently quantified in terms of the square of the Brunt-Väisälä frequency, ω 2 B (rad 2 s −2 ): Correspondence to: D. A. Hooper (d.a.hooper@rl.ac.uk) where θ (K) is potential temperature [=T (1000/p) 2/7 ], T (K) is absolute temperature, p (hPa) is pressure, g (m s −2 ) is gravitational acceleration and z (m) is altitude.Negative values of ω 2 B imply that the atmosphere is statically (convectively) unstable, whereas positive values imply that a vertically displaced air parcel will oscillate at frequency ω B .As will be seen shortly, the value of ω 2 B is typically small throughout the troposphere, ∼1×10 −4 rad 2 s −2 , it increases sharply at the tropopause level, and is ∼4×10 −4 rad 2 s −2 throughout the more statically stable lower-stratosphere.Perturbations from these values are commonly seen in association with a variety of atmospheric phenomena, including fronts (e.g.Röttger, 1979;Larsen andRöttger, 1983, 1985), gravity waves (e.g.Kitchen and Shutts, 1990) and turbulence (e.g.Browning and Watkins, 1970).Moreover, the variations of ω 2 B as a function of altitude have an effect on both gravity wave propagation (e.g.Prichard et al., 1995) and turbulence generation (e.g.Fritts and Rastogi, 1985).Routinely available altitude profiles of ω 2 B will therefore have a variety of scientific applications.
Balloon-borne radiosondes measure both temperature and pressure, at vertical intervals of a few tens of metres, and so the data are an ideal source of high-resolution ω 2 B profiles.However, radiosondes are seldom launched more than twice daily, whereas the time scales over which ω 2 B changes occur can be considerably less than 12 hours (e.g.Kitchen and Shutts, 1990).Moreover, radiosondes can drift horizontally more than 100 km from their launch site by the time that they reach an altitude of 15 km.As will be shown in Sect.5, ω 2 B values can change significantly within horizontal scales of just a few tens of kilometres.Profiles of radiosonde-derived ω 2 B cannot therefore necessarily be assumed to be representative of the atmosphere immediately above the balloon launch site.
Mesosphere-Stratosphere-Troposphere (MST) radar observations made in a vertically pointing direction (henceforth referred to as vertical beam observations) are another potential source of ω 2 B profiles.The frequency-power spectra of the radial velocity fluctuations often show a peak at ω B and a D. A. Hooper et al.: Static stability retrieval sharp decrease in power at higher frequencies (Röttger, 1980;Fukao et al., 1981;Ecklund et al., 1985;Revathy et al., 1996).It is therefore possible to derive a profile of ω 2 B with an altitude resolution equal to that of the radar observations -typically 150 or 300 m.A limitation of this technique is that it is necessary to make vertical beam observations at intervals of less than 1 min; the Brunt-Väisälä period is typically 10 min in the troposphere and 5 min in the lower stratosphere.For MST radars which are operated in the Doppler beam swinging mode, the interval between vertical beam observations is typically greater than 1 min.However, it is possible to overcome this by making multiple vertical beam observations within each radar cycle.A more serious limitation of this technique is that a clear spectral peak might not be apparent under conditions of large horizontal wind speeds (Ecklund et al., 1985).Scheffler and Liu (1986) showed that this was consistent with the effects of Doppler shifting by the horizontal winds.The technique is therefore probably of limited use at mid-latitude locations where upper-tropospheric wind speeds are typically several tens of m s −1 .
The present paper presents an alternative method which makes use of the vertical beam radar return signal power, P (Green and Gage, 1980).The two principal mechanisms responsible for lower-VHF radar returns are Fresnel scatter (Gage et al., 1981) and turbulent backscatter (Tatarski, 1961;VanZandt et al., 1978).The theories for both predict a direct relationship between P and ω 2 B through consideration of the square of the mean vertical gradient of generalised potential refractive index, M 2 (Ottersten, 1969).Although the relationship between P and M 2 has been demonstrated in a number of radar-radiosonde comparisons (e.g.Röttger, 1979;Green and Gage, 1980;Larsen andRöttger, 1983, 1985;Tsuda et al., 1988;Low et al., 1998), until now it has not been exploited for the retrieval of ω 2 B profiles.The present technique relies on finding an empirical fit between values of radiosonde-derived ω 2 B and of a radar-derived factor.The layout of the paper is as follows.A brief description of the radar and radiosonde observations is given in Sect. 2. The theoretical background to the technique is described in Sect.3. The linear regression analysis used to determine the best fit between the radar and radiosonde data is described in Sect. 4. The strengths and limitations of the technique are discussed in Sect. 5. Conclusions are given in Sect.6.

Radar and radiosonde observations
The data considered in this study are taken from the 52 MHz MST radar located at Esrange (Chilson et al., 1999), the Swedish Space Corporation's rocket range in northern Sweden (67.9 • N, 21.1 • E).Comparisons are made with data from radiosondes launched from the same site during the four consecutive winter periods (November-April) from 1996/1997-1999/2000.The Esrange MST radar has a peak transmitter power of 72 kW, a maximum duty cycle of 5%, and an antenna composed of a 140 5-element Yagi aerials spread over an area of 45×45 m.During the period in question, observations were typically switched, in sequence, between a variety of formats, which made use of both the spaced antenna and the Doppler beam swinging modes and a number of range resolutions.In the present study, attention will be confined to the format which is common to all periods: using the spaced antenna mode and with a range resolution of 300 m.Observations were made over the altitude interval 1.0-15.7 km at 300-m intervals, giving a total of 51 range gates.The data acquisition period for each observation was slightly less than 1 min, and the interval between observations varied depending on the other formats used within the sequence.
The radiosondes were of type Loran-C and had a typical ascent rate of 5 m s −1 .Data were recorded at 10-s intervals, giving a typical vertical separation of 50 m between consecutive measurements of temperature, pressure, relative humidity, wind speed and wind direction.The data were smoothed with a 300-m running mean, in order to give them an altitude resolution comparable to that of the radar data.Values of ω 2 B were subsequently calculated by taking gradients of potential temperature over 300-m intervals.The data were then transferred to the same altitude grid used for the radar observations by linearly interpolating between adjacent values, as necessary.The radiosondes typically took one hour to reach an altitude of 15.7 km.The radar return signal powers are therefore averaged over one hour, starting from the time of the radiosonde launch.

Theoretical background
In general, both Fresnel scatter and turbulent scatter can be expected to contribute to the radar returns observed by a vertically directed MST radar beam.One or the other mechanism will tend to dominate for specific time-altitude regions.The theory of Fresnel scatter (Gage and Balsley, 1980;Gage et al., 1981Gage et al., , 1985) ) predicts the following relationship between lower-VHF radar return signal power for a vertically directed beam, P (arbitrary linear units), and the mean vertical gradient of generalised potential refractive index, M (Ottersten, 1969): where α is a radar efficiency factor, P t is the peak transmitter power, A e is the effective area of the antenna, z is the range resolution, λ is the radar wavelength, and F (λ) is a wavelength-dependent constant of proportionality which relates the magnitude of the λ/2 harmonic component of M over the altitude interval z.It is noted that this expression has a z dependence (Gage et al., 1985), correcting a z 2 term used in the original formulation (Gage et al., 1981).Nevertheless, since Eq.(2) will only be considered in the context of a given radar, which operates at a fixed range resolution and with a fixed peak transmitter power, the relationship can be simplified to P ∝M 2 /z 2 .
For the case of turbulent scatter (Tatarski, 1961;VanZandt et al., 1978;Gage and Balsley, 1980) the following relationship is expected: where α is a ratio of eddy diffusion coefficients for potential refractive index and heat (this varies slightly with static stability but can be taken to be unity) and L 0 is the outer scale length of the turbulence spectrum.Again, this simplifies to P ∝M 2 /z 2 if it can be assumed that L 0 is a constant with respect to both time and altitude (VanZandt et al., 1978;Hooper and Thomas, 1998).This last assumption will be discussed in more detail once the retrieval algorithm has been outlined.
M is often conveniently approximated by the dry term, M D : although the full term additionally depends on both the specific humidity, q (kg kg −1 ), and its vertical gradient: It can be seen that the dry approximation is valid when the moduli of the second and third terms in the square brackets of Eq. ( 5) are significantly less than 1.The common assumption that this is the case above the first few kilometres of the atmosphere will need to be examined in more detail shortly.
For the time being the humidity contributions will be ignored entirely.
It will be recognised that the p/T term in Eq. ( 4) is proportional to density, which can be approximated as ρ 0 exp(−z/H ), where ρ 0 is the density at mean sea level and H (m) is the mean scale height across the considered altitude range (1.0-15.7 km).The local value of H is given by RT /g, where R is the gas constant for dry air.For the data considered in present investigation, the value of H is typically around 8 km at an altitude of 1 km and decreases to around 6 km at the tropopause level; it remains approximately constant at this value between the tropopause and an altitude of 15.7 km.A mean value of 6.71 km, calculated over all available data points, is adopted for H . Making the appropriate substitutions into P ∝M 2 /z 2 leads to the following simplified expression: which can be rearranged in order to define a radar factor r 2 B : where the altitude above mean sea level, z km , and the scale height, H km , are both given in units of kilometres, such that the following linear relationship is expected: The values of constants ω 2 0 and g 0 must be derived empirically from the best fit between radar-derived values of r 2 B and radiosonde-derived values of ω 2 B .This gives a retrieval algorithm which relies solely on parameters available from the radar, i.e. on z km and P .It is noted that the technique described here is essentially the same as that used by Gage and Green (1982) to determine temperature profiles, i.e. they integrated the profile of retrieved potential temperature gradient with respect to altitude.
The validity of the technique will be limited to those altitudes at which the humidity contributions to M can be ignored.The radar-derived values of ω 2 B will, in fact, be overestimated by a factor |M/M D |, i.e. by the modulus of the expression in the square brackets of Eq. ( 5).Since the necessary information is available from the radiosonde measurements, humidity-corrected values of r 2 B and of radar-derived ω 2 B , i.e. those multiplied by a factor |M D /M|, will also be considered in the following sections.It is noted that techniques exist for the retrieval of humidity information from radar return parameters (e.g.Stankov et al., 2003;Furumoto et al., 2003).However, these are beyond the scope of the current paper.
A final word should now be given to the necessary assumption that L 0 is a constant, with respect to both time and altitude, which is required to make the above derivation consistent with the theory of turbulent scatter.The equation for radar return signal power appears in the literature in a variety of alternative forms to that shown in Eq. ( 3).This is largely because the factor L 0 itself depends (inversely) on ω 2 B , leading several authors (e.g.Gage et al., 1980;Tsuda et al., 1988;Hocking and Mu, 1997;Hooper and Thomas, 1998) to suggest that P should be proportional to ω 2 B rather than to ω 4 B , as shown in Eq. ( 6).Moreover, Eq. (3) should contain an additional term describing the proportion of the radar observation volume containing turbulence, which again depends (inversely) on ω 2 B (VanZandt et al., 1978).Nevertheless, Hooper and Thomas (1998) demonstrated that profiles of vertical beam radar return signal power, at lowerstratospheric altitudes, were qualitatively related to those of ω 2 B regardless whether Fresnel scatter or turbulent scatter was the dominant mechanism.No attempt will be made to distinguish between the two types of radar returns in the current work.The proposed model, which is consistent with the relationship between radar return signal power (for both vertical and off-vertical beam observations) and M 2 demonstrated in a large number of radar-radiosonde comparisons (e.g., Röttger, 1979;Larsen andRöttger, 1983, 1985;Tsuda et al., 1988;Hocking and Mu, 1997;Hooper and Thomas, 1998;Low et al., 1998), will be justified by the results.

Results
The relationship between radiosonde-derived values of ω 2  to tropospheric measurements, and those with larger values of both, which correspond to lower-stratospheric measurements.Despite the fact that there is clearly a linear relationship between the two sets of values, there are also considerable deviations from this and the correlation coefficient is only 0.52.
Many of the largest deviations from a linear relationship, for the tropospheric measurements, can be accounted for by the fact that the humidity contributions to M have been ignored in the retrieval model.Figure 1 (right panel) shows the same data after humidity-correction has been applied to the r 2 B values.This correction has two main effects.Firstly, it shifts the cluster of data points corresponding to tropospheric values slightly to the left.Secondly, it shifts many of the outlying data points, with small values of ω 2 B and large values of of r 2 B , considerably to the left.The correlation coefficient increases from 0.52 to 0.68.It should be noted that the best fit between the data sets should be estimated using the humidity-corrected values of r 2 B .If the uncorrected values are used, the radar-derived values of ω 2 B will be overestimated at tropospheric altitudes for which the humidity contributions to M can be ignored.
The best fit is estimated following the method of Hocking et al. (2001) for comparing data sets measured by different techniques.For convenience the general symbols x and y will be used to refer to the humidity-corrected radar-derived values of r 2 B and the radiosonde-derived values of ω 2 B , respectively.The lines with gradients g x and g y , shown in Fig. 1 (right panel), represent, respectively, the least-squares best fits from the regression of y on x and vice versa.These correspond, respectively, to the assumptions that all of the variability is associated with the y values, and that there are no errors associated with the x values, and vice versa.In addition to the fact that there are measurement errors associated with both the radar and the radiosonde data, each instrument is measuring different parameters, and each parameter relates to different temporal and spatial scales of atmospheric variability (Kitchen, 1989).Moreover, the horizontal separation between the radar and radiosonde measurements increases with increasing altitude.As will be discussed in the next section, it appears that this last factor is responsible for the largest degree of variability between the two data sets.Although neither of the lines shown in Fig. 1 (right panel) represents the most appropriate fit between the two data sets, each one represents an extreme of the possible fits; the required fit lies somewhere between the two.
In order to find these extreme fits, it is first necessary to discard a small number of outlying data points which would otherwise have a disproportionately significant influence.
Starting with all 11 050 data points, the least-squares fitting routine is performed recursively, removing data points lying more than 3 standard deviations from the fit at each step, until no more data points are removed.This is performed separately for the regression of y on x, which leaves 10 699 data points after 10 steps, and for the regression of x on y, which leaves 10 621 data points after 9 steps.The lines shown in Fig. 1  10 527 remaining data points which are common to both removal procedures.Discarding less than 5% of the original data points has the effect of increasing the correlation coefficient from 0.68 to 0.82.
Standard deviations, σ x and σ y , which relate to the radar and radiosonde measurements, respectively, contain information about the intrinsic measurements errors, about the natural spatial and temporal variability of the atmosphere, and about the design of the experiment being used to perform the measurements.The values of σ x and σ y are uniquely interrelated with that of the best fit gradient g 0 , as shown in Fig. 2; see Hocking et al. (2001) for full details.The possible values of σ x and σ y are therefore constrained by the fact that the value of g 0 lies somewhere in the range g x (5.97×10 −8 ) to g y (8.94×10 −8 ).
Since it is thought that the horizontal separation between the radar and radiosonde measurements is responsible for the largest deviations from a linear relationship (this will be justified in the following section), the value of g 0 is chosen such that the variability is equally distributed between the radar and radiosonde measurements, i.e. to the condition that σ y =g 0 σ x .It can be seen from Fig. 2 that this results in g 0 =7.30×10 −8 , σ x =9.75×10 2 arbitrary radar units, and σ y =7.12×10 −5 rad 2 s 2 .The values of σ x and σ y set upper limits on the measurement errors associated with the radar and radiosonde data, respectively.The intercept of the best fit line with the y-axis (ω 2 0 =−5.44×10−6 rad 2 s −2 ) is chosen such that the best fit line and the lines with gradients g x and g y all intersect at one point.

Figure 3 compares profiles of ω 2
B derived from radiosonde data (thick grey lines) and radar data (thick black lines) for 5 individual cases.The profiles are clearly quantitatively and qualitatively well matched.The only significant discrepancies occur in the lower-troposphere.As can be seen from the corresponding profiles of humidity-corrected radar-derived values (thin black lines), these can be attributed to the fact that humidity contributions to M have been ignored in the retrieval model.This is clearly the case at altitudes below 4 km in the fourth and fifth panels from the left in Fig. 3. On the other hand, the portions of the radar-derived profiles with values of less than 1×10 −4 rad 2 s −2 , such as those seen between the altitudes of 4 and 8 km in the middle panel of Fig. 3, can be assumed to be representative, even if slightly quantitatively inaccurate.
It is noted that the use of humidity-correction, applied to radar-derived values of ω 2 B , does not always lead to the improvements in fit to the radiosonde-derived values as suggested by the cases shown in Fig. 3. Similar results have been reported in a number of other studies and the discrepancies are known to be particularly large under conditions of precipitation (e.g.Chu and Lin, 1994;Vaughan and Worthington, 2000).
Not all of the discrepancies between radar-and radiosonde-derived values of ω 2 B are caused by limitations of the retrieval technique.Some are simply caused by the fact that the radar and radiosondes are sampling different regions of the atmosphere.For the example shown in Fig. 5, the radiosonde is over 46 km away from the radar by the time that it reaches an altitude of 9 km, over 77 km away at 13 km altitude, and over 90 km away at 16 km altitude.For the region between 9 and 13 km, the profiles of ω 2 B (left panel) are poorly matched both quantitatively and qualitatively.
The fact that the perturbations of radiosonde-derived ω 2 B , between 9 and 13 km altitude, are anti-correlated with those of the radiosonde's ascent rate (right panel), is a clear indication that the radiosonde is passing through a region of significant mountain wave activity.Mountain wave induced perturbations of ω 2 B are not generally apparent at tropospheric altitudes, although the perturbations of ascent rate suggest that the radiosonde is experiencing mountain wave activity throughout the depth of the troposphere.Kitchen and Shutts (1990) showed that the amplitude of such ascent rate perturbations typically increases exponentially with increasing altitude.They attributed sharp decreases in amplitude, such as that seen at 13 km altitude in the present example, to the radiosonde passing out of a horizontally localised region of wave activity.They found that significant changes in temperature perturbation amplitude, at a given altitude, could occur over horizontal distances of just 10 km.In the present example, the wave activity reappears, and persists, when the radiosonde is above 16 km altitude (not shown) and over 90 km away from the radar.This example highlights the fact that radiosonde measurements are not necessarily representative of the atmosphere directly above the launch site.

Conclusions
This paper has demonstrated that the well-known relationship between the radar return signal power (for a vertically directed beam), P , and the mean vertical gradient of generalised potential refractive index, M, can be exploited in order to retrieve profiles of the square of the Brunt-Väisälä frequency, ω 2 B , from lower-VHF radar observations of the atmosphere.The technique relies on fitting a radar-derived factor, which principally depends on the square root of P , to radiosonde-derived values of ω 2 B .The technique is only valid for those altitudes at which humidity contributions to M can reasonably be ignored.In general, this limits the applicability to lower-stratospheric altitudes, although under specific conditions it can be used over sections of the troposphere.A strength of the technique is that the derived profiles represent the nature of the atmosphere directly above the radar.Since radiosondes drift horizontally as they rise through the atmosphere, the data are not necessarily representative of the atmosphere above the launch site.This is shown to be particularly significant when the radiosonde experiences mountain wave activity, which can change considerably over horizontal distances of just 10 km.

BFig. 1 .
Fig. 1.The relationship between radiosonde-derived values of ω 2B and radar-derived values of r 2 B .The plot area is divided into a 100 by 100 grid, with divisions of 0.1×10 −4 rad 2 s −2 along the y-axis and of 0.1×10 3 arbitrary r 2 B units along the x-axis; the grey scale represents the number of data points falling within each cell.The values of r 2 B shown in the right-hand panel have been humidity-corrected.The lines labelled g x and g y represent the best fits from the regression of y on x and vice versa.

Fig. 5 .
Fig. 5.An example of large discrepancies between the radiosondederived (thick grey line) and radar-derived (thick black line) values of ω 2 B (left panel) corresponding to 08:18 UT on 1 November 1999.The radiosonde's ascent rate is shown in the right panel.