26 Feb 2021
26 Feb 2021
Fourdimensional mesospheric and lower thermospheric wind ﬁelds using Gaussian process regression on multistatic specular meteor radar observations
 ^{1}Haystack Observatory, Massachusetts Institute of Technology, USA
 ^{2}Leibniz Institute of Atmospheric Physics at the University of Rostock, Germany
 ^{3}UiT Arctic University of Norway, Norway
 ^{1}Haystack Observatory, Massachusetts Institute of Technology, USA
 ^{2}Leibniz Institute of Atmospheric Physics at the University of Rostock, Germany
 ^{3}UiT Arctic University of Norway, Norway
Abstract. Mesoscale dynamics in the mesosphere and lower thermosphere (MLT) region have been difficult to study from either ground or satellitebased observations. For understanding of atmospheric coupling processes, important spatial scales at these altitudes range between tens to hundreds of kilometers in the horizontal plane. To date, this scale size is challenging observationally, and so structures are usually parameterized in global circulation models. The advent of multistatic specular meteor radar networks allows exploration of MLT mesocale dynamics on these scales using an increased number of detections and a diversity of viewing angles inherent to multistatic networks. In this work, we introduce a four dimensional wind field inversion method that makes use of Gaussian process regression (GPR), a nonparametric and Bayesian approach. The method takes measured projected wind velocities and prior distributions of the wind velocity as a function of space and time, specified by the user or estimated from the data, and produces posterior distributions for the wind velocity. Computation of the predictive posterior distribution is performed on sampled points of interest and is not necessarily regularly sampled. The main benefits of the GPR method include this nongridded sampling, the builtin statistical uncertainty estimates, and the ability to horizontallyresolve winds on relatively small scales. The performance of the GPR implementation has been evaluated on Monte Carlo simulations with known distributions using the same spatial and temporal sampling as one day of real meteor measurements. Based on the simulation results we find that the GPR implementation is robust, providing wind fields that are statistically unbiased and with statistical variances that depend on the geometry and are proportional to the prior velocity variances. A conservative and fast approach can be straightforwardly implemented by employing overestimated prior variances and distances, while a more robust but computationally intensive approach can be implemented by employing training and fitting of model parameters. The latter GPR approach has been applied to a 24hour data set and shown to compare well to previously used homogeneous and gradient methods. Small scale features have reasonably low statistical uncertainties, implying geophysical wind field horizontal structures as low as 20–50 km. We suggest that this GPR approach forms a suitable method for MLT regional and weather studies.
 Preprint
(19718 KB) 
Supplement
(9531 KB)  BibTeX
 EndNote
Ryan Volz et al.
Status: final response (author comments only)

RC1: 'Comment on amt202140', Anonymous Referee #1, 15 Mar 2021
This is an important, clearlywritten article which describes the
application of Gaussian process regression (GPR) to MLT meteor radar
winds data acquired with multistatic MIMO meteor networks which
produce angle of departure and angle of arrival observables in
addition to Doppler shifts. Thanks to wavevector diversity, data from
the network can be used to characterize the threedimensional wind
field without the ambiguity inherent in individual monostatic systems
which provide no information about horizontal vorticity. The GPR
method is tested here as a framework for resolving finescale
structure in the flow field and for performing error propagation. The
GPR formalism is developed and tested using synthetic and actual
data. The results are compared with those from the more conventional
homogeneous and gradient methods.Ovearll, the significance of the paper is excellent in view of the
importance of the MLT winds (which are crucial for understanding
ionosphere/atmosphere coupling) but difficult to measure using other
groundbased methods. Also noteworthy is the novelty of the
multistatic method which is revitalizing the meteorradar field. The
potential of the method can be realized through incisive analysis such
as is provided by GPR. The scientific method is also excellent in view
of the appropriateness of the methodology and the clarity of the
fomalism and explanations. The presentation quality is good and could
be improved by clarifying the meaning of some of the figures.The paper is acceptible for publication but could benefit from
addressing a few minor questions:1. Throughout the paper, the GPR method is described as being
nonparametric. However, the prior model for the winds is necessarily
parametrized by the means and variances of the wind components
(theta). In what sense is this method nonparametric?2. The paper distinguishes the GPR method from Tikhonov regularization
which is viewed as a competing method. However, 2nd order Tikhonov
regularization can be interpreted as the adoption of a Gaussian prior
for the state vector. How is this fundamentally different from the
method presented here?3. The paper describes the assumption of a Gaussian random process for
the winds as "convenient" because of the tractable computations that
result. Could not the authors provide a more satisfying rationale by
considering whether the central limit theorem applies to the MLT wind
components?4. The authors of the paper note the absence of other kinds of MLT
wind measurements that could be used to validate the wind estimates
produced by the GPR method. Have the authors considered the use of
generalized cross validation?5. The authors note that the region of validity and the spatial
resolution of their method depends on the geometrical configuration of
the multistatic meteor array. Have they considered developing a
geometrical dilution of precision (GDOP) estimator?6. Finally, interpreting figures 24 and 7 is very difficult in view
of the fact that color gradations are being used to represent multiple
quantities simultaneously (i.e. horizontal winds, vertical winds, and
data quality). The authors should attempt to clarify these figures.
AC1: 'Reply on RC1', Ryan Volz, 30 Apr 2021
Thanks for your thoughtful comments. We will address each of them in our formal reply and revised version, but for now we would like to comment on the following selected points:
1. (GPR as "nonparametric") We describe GPR as nonparametric in the same sense as it is used in the referenced Rasmussen and Williams (2006) book: an estimation method which does not compress the training data into a finitedimensional parameter vector, in contrast to parametric methods like linear regression. The mean and covariance function parameters are usually called hyperparameters to emphasize that they are parameters of a nonparametric model. We also find this terminology confusing, so we will seek to clarify it in the next revision.
2. (Comparison to Tikhonov regularization) The methods are indeed related, and there is a good discussion of this topic throughout Rasmussen and Williams (2006), particularly in Chapter 6. The most direct difference is that Tikhonov regularization would best relate to GPR with a squared exponential covariance, whereas we have employed a Matern5/2 covariance. That detail strikes at the heart of the difference between the two methods: it can often be more natural to express prior knowledge in terms of a GPR covariance than as a likelihood penalization. Practically, the difference also comes down to how it is more natural to perform nongridded estimation and analyze uncertainty with GPR compared to regularization approaches. In many respects, the approaches are two sides of the same coin, which is why we see value in future intercomparison that can help refine both approaches.
5. (Geometrical dilution of precision estimator) Perhaps it is clearer to say that the measurement geometry controls the wind estimate uncertainty, which we can calculate through the posterior covariance, and naturally there are (location, wind direction) pairs that have higher uncertainty. As far as we understand it, we can quantify the GDOP as a function of location in this case by taking the root mean square of the wind component uncertainties. If this doesn't describe what you were pointing toward, then we'd love to hear more details!
6. (Figure clarity) We very much would like these figures to be both expressive and interpretable, and we recognize that this is a difficult task given the amount of information that they attempt to include. We strived to use color scales for the different elements that are distinct enough to be identifiably separate and thought we had achieved a good balance. Given that, we're uncertain about how exactly to make satisfying improvements, although we can try. Do you have specific detailed critiques that we could implement (e.g. the vertical wind colors on the streamlines are too hard to distinguish from the background colors, or the information conveyed with them is not worth the visual clutter)?

AC7: 'Further Reply on RC1', Ryan Volz, 16 Sep 2021
Thanks again for your thoughtful comments. Here are our comments on the remaining points:
3. (Justification for assuming Gaussian random processes) This question prompted us to think more about this assumption, so we thank you for that. Assuming normality imposes the minimal prior information about the wind processes within a secondorder statistical framework since the Gaussian distribution has maximum entropy for a known mean and covariance. We will add this explicit justification to the revised manuscript.
4. (Generalized crossvalidation) Generalized crossvalidation is indeed one technique that could be used to further assess the validity of the wind estimates using just the data we already have. There is a good discussion of crossvalidation in the context of GPR in Rasmussen and Williams (2006), section 5.4. We think this is a good area for future work (and will additionally note it in the revised manuscript), and the result would be a better understanding of the covariance hyperparameter values and the merit of different forms for the prior wind distributions. We leave it to future work because it addresses the topic of improving the model specification, and that area is large enough to merit its own paper(s) beyond the introduction of the technique that we provide here. Although such work would speak to validation, it still remains that comparison to a completely independent technique would go even further to give confidence in the GPR method.
6. (Figure clarity) Considering this and other reviewer and community comments, we will be simplifying the figures somewhat in the revised manuscript. Notably, we will be removing the vertical wind coloring from the horizontal wind streamline maps. We decided that the focus of these figures should be the horizontal winds, and including the additional color axis invites confusion and diverts focus. Additionally, we will be simplifying the bias panels in Figures 3 and 4 to remove the green shading and change the contour lines to depict the predictive variance values instead of dB improvement. This makes the comparison to the error variance more direct, while still giving the viewer the necessary information to focus on the more pertinent regions of the bias values.

AC1: 'Reply on RC1', Ryan Volz, 30 Apr 2021

CC1: 'Comment on amt202140', Chris Meek, 27 Mar 2021
The comment was uploaded in the form of a supplement: https://amt.copernicus.org/preprints/amt202140/amt202140CC1supplement.pdf

AC2: 'Reply on CC1', Ryan Volz, 01 May 2021
We appreciate your concern. We believe there is a misunderstanding, since our Dopplerderived lineofsight velocity measurements are no different than what have been used for decades to estimate winds from specular meteor scatter, including the Hocking et al. (2001) and Holdsworth et al. (2004) papers cited in the first paragraph of section 2. At least for an ideal meteor trail, we agree with the motion and measurement depicted in your figure. The trail moves from (A) to (B) with a given velocity, and we measure the projection of this velocity onto the line of sight which moves the reflection point from (A) to (C). Where we disagree is that this "creat[es] an apparent vertical motion" in the first case or "horizontal motion" in the second case. Yes, the projection has a vertical (horizontal) component, but that does not mean that we conclude from those measurements that the wind velocity has a vertical (horizontal) component. From a single measurement, we can't know what the full velocity vector is. But by combining nearby measurements with assumptions about the smoothness of the wind field and quantification of the measurement uncertainty, we can estimate the complete wind velocity vector. The measurement geometries provided by the multistatic configuration allow for such "overlapping" measurements of different meteors with a diversity of projection directions, so this works well in many cases. More importantly, we also know when we don't have sufficient information to resolve the full wind vector, and the estimate uncertainties reflect that. Indeed, the uncertainties for the vertical wind component are generally higher than the horizontal components relative to the mean because we observe relatively few meteor trails with a large vertical line of sight projection (trails close to horizontal). We hope this addresses your concern and clarifies the technique.

CC2: 'Reply on AC2', Chris Meek, 03 May 2021
Thank you; that was a good response to my comment and recognizes the problem,
but doesn't solve it. Meteor wind papers all seem to accept the handeddown wisdom
that the "echo" motion is equal to the lineofsight (LOS) component of
the wind. But based on this assumption, and that all echoes are seen above the
horizon, then a vertical wind must automatically exist, no matter what the trail
orientation (rotation around the LOS). The early papers were analysing for horizontal
wind, so the actual scatter mechanics were not as important.
Because I only have experience with a monostatic system some of my followjng comments
probably apply only approximately to a bistatic system. But hope that
these considerations are completely mitgated by it is not enough to convince.Because the actual expected wind vertical component are in cm/s rather
than m/s, 2D wind (horizontal) values are reasonable and agree well
with other measurements, and are correct under the assumption of
zero vertical velocity. If 3D fits are done, extreme vertical
wind values are often found  I have found nonspurious 2030 m/s. At the
time I had no explanation.
If meteor echo distribution were uniform in azimuth, the problem of vertical
wind artifacts could be mitigated to some extentf  but it is not
uniform. Statistically the azimuthal direction rotates during the day
(at least at Eureka). So a multistatic system unless maybe very large spacing
doesn't change this[ I tried once to fix this lopsided echo distribution by dividing into
octants and equalizing the weight given to each in the wind fit  but there were always
one or two octants with very few, or no, meteors ]
Another related comment (while I have the podium) is about component wind errors/perturbations.
I don't know if this applies to the current paper.
A standard least squares fit allows the formal calculation of component errors. But constant
wind models with uniform echo distribution and an artficial perturbation/error added
in one component show that the variation "bleeds" into the other component. That's because of the
radial nature of the sampling. A perturbation in zonal wind, say, causes a perturbation in
radial speed component, which appears as a perturbation in the meridional component.
I don't see a solution for these problems, and they must have some effect on results.
particularly as more detailed and complicated analyses are used. Caveats should be stated
(or assumptions, like zero vertical velocity  which I think is the hidden, but understood,
assumption in horizontal wind analysis).
Trail orientation can potentially be estimated from echo polarization  but I don't
think this helps here  though might be interesting for meteor studies.
Finally, I think confirmation of the results requires a realistic model, including a variety of
trail orientations and a sporadic meteor model (e.g. Margaret CamplbellBrown's)
 Unfortunately this is not a simple process for bistatic.
cem

AC3: 'Reply on CC2', Ryan Volz, 06 May 2021
There is a lot we agree about regarding the difficulties and pitfalls of the meteor wind measurement, but perhaps we're not quite speaking the same language. We'll try to clarify a few relevant points:
a) It is suitable to apply intuition from the monostatic case to the bistatic case, with minor adjustment to a notional equivalent monostatic system. For each meteor scattering, the bistatic scattering is equivalent to a monostatic geometry where the middle point of the bistatic link gives the virtual monostatic radar location, and the loci of zero velocity is an ellipse with foci at the tx and rx instead of a circle. The effective Bragg wavelength will depend on the geometry and will be greater than or equal to half the radar wavelength, and the "lineofsight" direction is given by the difference between the scattering and incident wave vectors. The particulars of the bistatic case are discussed in more detail in Stober and Chau (2015) and Chau et al. (2019).
Stober, G., & Chau, J. L. (2015). A multistatic and multifrequency novel approach for specular meteor radars to improve wind measurements in the MLT region. Radio Science, 50(5), 431–442. https://doi.org/10.1002/2014RS005591
Chau, J. L., Urco, J. M., Vierinen, J. P., Volz, R. A., Clahsen, M., Pfeffer, N., & Trautner, J. (2019). Novel specular meteor radar systems using coherent MIMO techniques to study the mesosphere and lower thermosphere. Atmospheric Measurement Techniques, 12(4), 2113–2127. https://doi.org/10.5194/amt1221132019
b) Accounting for Earth curvature is an important consideration when working with measurements over an area of this scale. This shows up in differences between what would be seen as "vertical" or "horizontal" motion in the EastNorthUp reference frame at the radar receiver versus the local vertical and horizontal directions in the ENU frame at the meteor location. There is a discussion of this and the coordinate conversion procedure that we use in Stober et al. (2018). Whenever we refer to vertical or horizontal winds, we mean with respect to the Earth's surface underneath that location taking into full account these coordinate conversions. Perhaps that addresses some of the concerns here.
Stober, G., Chau, J. L., Vierinen, J., Jacobi, C., & Wilhelm, S. (2018). Retrieving horizontally resolved wind fields using multistatic meteor radar observations. Atmospheric Measurement Techniques, 11(8), 4891–4907. https://doi.org/10.5194/amt1148912018
c) The LoS Doppler measurement indeed contains contributions of vertical velocity, and most works based on monostatic systems have ignored the vertical velocities because of sampling issues, inhomogeneous horizontal winds, etc. The assumption has been that horizontal winds are much larger than vertical winds, and the mean vertical wind is very small. The GPR technique is agnostic to the assumptions the user wants to impose about the horizontal winds, and one could force a zero vertical wind component or scale the variance of the vertical component relative to the horizontal to suitably affect the final estimated winds. There are important estimation decisions here that welcome further consideration, but we consider that to be beyond the scope for the introduction of the GPR technique as contained in this paper.
d) To that point, the Ph.D. student of one of the coauthors (JLC) is about to submit a paper dealing with the vertical velocity estimates using a virtual meteor radar system (i.e., geometry and realistic sampling) on a high resolution atmospheric model (Upper Atmosphere  ICON). Namely, the measured line of sight velocities are replaced by projected velocities using UAICON winds at each meteor detection point. One of the main conclusions is that a variability of ± 12 m/s is an apparent vertical wind variability due to horizontal wind inhomogeneities. So we agree, mean vertical wind estimates in monostatic and bistatic systems obtained over the whole area (a few hundred of kms diameter) present a large artificial variability.
e) In this work, our focus is in the wind estimates with much smaller horizontal scales than traditional meteor systems, taking advantage of more samples and different viewing angles (equivalent as having different monostatic systems with relatively close separations). That is contrary to previous works: we don’t do a global 3D wind fit on the illuminated area, but instead we effectively do it locally in smaller areas (controlled by the covariance length scale). We show that the reliable area of vertical velocities is narrower than the area for horizontal velocity estimates. These areas depend on the sampling geometry that in turns depend not only in the system geometry but also, as you mentioned, on the meteor population.
f) Although we don't highlight it, the GPR technique also provides the twocomponent covariance values as final estimation outputs in addition to the singlecomponent posterior variances. These covariance entries allow one to quantify exactly one of the issues you mention: how perturbation/error in one wind component bleeds into the estimate of another wind component. This doesn't solve the problem of correlated errors, but it does let us be aware of it, quantify it, and potentially optimize meteor radar network geometries to minimize it. This is an important topic that begs for further investigation, and we think that the GPR outputs/analysis can help with that.

AC3: 'Reply on CC2', Ryan Volz, 06 May 2021

CC2: 'Reply on AC2', Chris Meek, 03 May 2021

AC2: 'Reply on CC1', Ryan Volz, 01 May 2021

CC3: 'Comment on amt202140', Gunter Stober, 18 Jun 2021
This is an interesting approach for multistatic meteor radar networks. Some part of the retrievals might deserve some additional information about the inversion method:
1. How does the method differ from a classical 2DVAR data assimilation approach applying a standard cost function? The inversion is computed layerbylayer using a 'large vertical averaging'? However, the atmosphere shows very often strong vertical shears. Are these shears considered in the cost function or error estimates or both?
2. The algorithm makes use of the WGS84 geometry. It might be appropriate to cite the original paper and references of the implemented algorithms.
3. Retrieving vertical winds is usually challenging due to exponential instability growth for parameters with low measurement response. This often limits the vertical resolution for radiometers. Do they check the solutions for such instabilities? The values that can be found in the examples are rather high.

AC4: 'Reply on CC3', Ryan Volz, 24 Jun 2021
Thanks for the comment!
1. (Comparison to classical 2DVAR data assimilation) Can you point to a reference to which we might compare? In general we expect similarities with many statistical estimation approaches, with the primary differences arising from the specific formalism used and how that allows one to express the prior knowledge and constraints on the estimate. We think the strengths of the GPR method are that specifying a covariance function allows one to easily apply physical intuition, that it works in 4 dimensions at once and does not impose a gridding of the measurements or estimates, and that it is natural to carry through and compute uncertainties on the estimates.
In particular concerning vertical averaging and accomodating strong vertical shears, the natural way to express that in the GPR framework is to enforce a small value for the covariance length scale in the vertical dimension. The examples in this paper use a fitted value for the vertical length scale of 3 km, which arises both from what the density of meteor measurements will support and also from the observed vertical correlations particular to the winds seen with these data. So we can say that the averaging is approximately over a vertical region of 3 km and that any sharper vertical shears will be smoothed over (i.e. the effective vertical resolution is 3 km). It is certainly possible to manually set the covariance parameters to better investigate vertical shears in particular if that is what the user sets out to do. The GPR framework provides great flexibility, and we certainly don't think that we have already arrived at a form for the covariance functions that achieves best case performance in general and especially for specific analysis objectives.
2. (WGS84 geometry) There are two relevant applications of applying WGS84 geometry related to this work: in geolocating the meteors and calculating the corresponding Bragg vector in the meteor local ENU coordinate system, and in estimating the winds from the meteor data. This paper only lightly touches on the first case in assuming that the meteor radar data has already been processed into meteor observations complete with a Doppler shift, geodetic latitude and longitude, and Bragg vector. The procedure to get to that point is exactly the same as in the works cited in Section 4, e.g. Chau and Clahsen (2019), although full detail can be found in Clahsen (2018) [full citation below], with the meteor geolocation procedure also described in Stober et al. (2018). We will clarify this in the next revision.
For the second case in applying WGS84 geometry within the wind estimation, we describe our procedure of using a local azimuthal equidistant projection in Section 2. This application is new to this work, although the azimuthal equidistant projection using the WGS84 sphereoid is wellknown. We use a projection only for computational efficiency, as it is much easier and faster to work in a Euclidean coordinate system when performing the estimation than to use the full geodetic computations, although it is possible to do the latter.
3. (Vertical winds) We agree on the challenges of estimating vertical winds and the care that must be taken in doing so. This procedure takes into account the projection effects that result in the vertical wind component usually not being as wellcharacterized as the horizontal components, and this relative uncertainty is quantified in the posterior estimated variances. This is best demonstrated in the simulation results and in particular Figure 4. In general, the user of this technique can impose their prior belief in generally small vertical wind values by appropriately setting a small value for the vertical wind covariance amplitude or by checking that fitted values correspond to that expectation.
There are a few locations in the examples where the estimated vertical wind exceeds a magnitude of 10 m/s. We don't make a note of it because the estimated confidence interval of the vertical wind at those points (and indeed, all over) is large enough to allow for a true value that is in line with what has been observed in the MLT region by different instruments and groups, e.g. Vierinen et al. (2013), Lu et al. (2017), Chau et al. (2020). These examples show MLT vertical winds on the scale of a few meters per second at minute to hour timescales.
Additional ReferencesChau, J. L., Urco, J. M., Avsarkisov, V., Vierinen, J. P., Latteck, R., Hall, C. M., & Tsutsumi, M. (2020). FourDimensional Quantification of KelvinHelmholtz Instabilities in the Polar Summer Mesosphere Using Volumetric Radar Imaging. Geophysical Research Letters, 47(1), e2019GL086081. https://doi.org/10.1029/2019GL086081
Clahsen, M. (2018, January 31). Error analysis of wind estimates in specular meteor radar system (M.S. Thesis). University of Rostock, Rostock, Germany.
Lu, X., Chu, X., Li, H., Chen, C., Smith, J. A., & Vadas, S. L. (2017). Statistical characterization of hightomedium frequency mesoscale gravity waves by lidarmeasured vertical winds and temperatures in the MLT. Journal of Atmospheric and SolarTerrestrial Physics, 162, 3–15. https://doi.org/10.1016/j.jastp.2016.10.009
Vierinen, J., Kero, A., & Rietveld, M. T. (2013). High latitude artificial periodic irregularity observations with the upgraded EISCAT heating facility. Journal of Atmospheric and SolarTerrestrial Physics, 105–106, 253–261. https://doi.org/10.1016/j.jastp.2013.08.012

CC4: 'Reply on AC4', Gunter Stober, 25 Jun 2021
Dear Ryan,
thanks for clarifying the 2D nature of the presented retrievals. From the reply I take that the winds are computed within a 2D layer with an irregular horizontal grid using a 3 km vertical average. This setup sounds reasonable. The terms 2DVAR or 3DVAR or higher are typically used for data assimilation e.g., Gelaro et al., 2017 for MERRA2 (3DVAR) and Eckerman et al., 2019 for NAVGEMHA (4DVAR) and references therein. It might be worth describing the algorithm concerning this aspect with some more details.
WGS84 geometry implementation:
The reply about the WGS84 geometry implementation somehow confused me. Chau and Clahsen (2019) cite Stober et al., 2018 as a source of the WGS84 geometry. I am even more confused as one of the coauthors of this paper (Jorge L. Chau) confirmed this in a public comment a few weeks ago (see link: https://doi.org/10.5194/amt2021124CC2).
I might also want to recall that we had a good meeting in Kühlungsborn in March 2017. The attendees were Ryan Volz, Philip J. Erickson, Juha Vierinen, Jorge L. Chau and Gunter Stober. We discussed the multistatic projects and presented the status and help the MIT colleagues to catch up with the multistatic developments that were done so far until 2017. The patent is mainly based on the details presented in Stober et al., 2018 including the WGS84 geometry, which was already in parts implemented in Stober and Chau (2015) , but only to solve the forward scatter geometry and the altitude of the meteors without angular correction. Due to the embargo time between patent submission and the patent acceptance the submission to AMT had to be delayed. If the editor needs further evidence for the timeline, I can forward the documents in a private email to avoid public embarrassment.
I also tried to google the cited thesis of Clahsen (2018), but I received a match to another AMT paper, which cites Stober et al., 2018.
Vertical winds:
I am also confused about the reply on the vertical winds. They now justify the vertical winds as real by presenting additional citations. This is contradicting another time a public comment by Jorge L. Chau (see again link from above) related to the SIMONe systems. In the public comment, Jorge L. Chau clearly confirms that there is no claim that the winds are physically meaningful or presenting the geophysical truth. I suggest they add a paragraph discussing this obvious contradiction or they present a physical explanation within the laws of thermodynamics that these winds are realistic. Considering that the residual circulation causes a 100 K deviation from the thermal equilibrium at the mesosphere lower thermosphere (Becker, 2012, Smith, 2012) and is associated to mean vertical velocities of a few mm/s to cm/s. So it might be worth estimating the adiabatic cooling of the presented 10 m/s upwelling and estimating the cooling rate per day. If these winds are real, this will be a gamechanger of our understanding of the middle atmosphere dynamics at the MLT. Otherwise, they can also discuss the obvious controversy to other observations. I added some references, which show different vertical velocities. It would be good to discuss these things a bit more. It is the one point that doesn’t fit, that brings science forward.
I really enjoyed reading the paper and I am looking forward to future scientific results.
Gelaro, R., McCarty, W., Suárez, M. J., Todling, R., Molod, A., Takacs, L., Randles, C. A., Darmenov, A., Bosilovich, M. G., Reichle, R., Wargan, K., Coy, L., Cullather, R., Draper, C., Akella, S., Buchard, V., Conaty, A., da Silva, A. M., Gu, W., Kim, G., Koster, R., Lucchesi, R., Merkova, D., Nielsen, J. E., Partyka, G., Pawson, S., Putman, W., Rienecker, M., Schubert, S. D., Sienkiewicz, M., & Zhao, B. (2017). The ModernEra Retrospective Analysis for Research and Applications, Version 2 (MERRA2), Journal of Climate, 30(14), 54195454.
Eckermann, S. D., Ma, J., Hoppel, K.W., Kuhl, D. D., Allen, D. R.,Doyle, J. A., Viner, K. C., Ruston, B. C., Baker, N. L., Swadley, D., Whitcomb, T. R., Reynolds, C. A., Xu, L., Kaifler, N., Kaifler,B., Reid, I. M., Murphy, D. J., and Love, P. T.: HighAltitude (0–100 km) Global Atmospheric Reanalysis System: Description and Application to the 2014 Austral Winter of the Deep Propagating Gravity Wave Experiment (DEEPWAVE), Mon. Weather Rev., 146, 2639–2666, https://doi.org/10.1175/MWRD170386.1, 2018.
Smith, A. K.: Global Dynamics of the MLT, Surv. Geophys., 33, 1177–1230, https://doi.org/10.1007/s1071201291969, 2012.
Becker, E.: Dynamical Control of the Middle Atmosphere, Space Sci. Rev., 168, 283–314, https://doi.org/10.1007/s1121401198415, 2012.
Vertical wind observations:
Most of them contradict the wind speeds cited in the manuscript (certainly not complete).
Chau, J. L., Stober, G., Hall, C. M., Tsutsumi, M., Laskar, F. I., and Hoffmann, P.: Polar mesospheric horizontal divergence and relative vorticity measurements using multiple specular meteor radars, Radio Science, 52, 811–828, ttps://doi.org/10.1002/2016RS006225, 2017.
Fritts, D., Hoppe, U.P., and Inhester, B.: A study of the vertical motion field near the highlatitude summer mesopause during MAC/SINE, J. Atmos. Terr. Phys., 52, 927–938, https://doi.org/10.1016/00219169(90)90025I, 1990.
Hoppe, U.P. and Fritts, D. C.: Highresolution measurements of vertical velocity with the European incoherent scatter VHF radar: 1. Motion field characteristics and measurement biases, J. Geophys. Res., 100, 16813–16825, https://doi.org/10.1029/95JD01466, 1995a.
Straub, C., Tschanz, B., Hocke, K., Kämpfer, N., and Smith, A. K.: Transport of mesospheric H_{2}O during and after the stratospheric sudden warming of January 2010: observation and simulation, Atmos. Chem. Phys., 12, 5413–5427, https://doi.org/10.5194/acp1254132012, 2012.
I. Laskar, J. L. Chau, J. P. StMaurice, G. Stober, C. M. Hall, M. Tsutsumi, J. Höffner und P. Hoffmann, Experimental evidence of Arctic summer mesospheric upwelling and its connection to cold summer mesopause, Geophys. Res. Lett., 44, 91519158, doi:10.1002/2017GL074759, 2017.
Gudadze, G. Stober und J. L. Chau, Can VHF radars at polar latitudes measure mean vertical winds in the presence of PMSE?, Atmos. Chem. Phys., 19, 44854497, doi:10.5194/acp1944852019, 2019.
Vincent, R. A., Kovalam, S., Murphy, D., J, Reid, I. M., & Younger, J. P. (2019). Trends and variability in vertical winds in the southern hemisphere summer polar mesosphere and lower thermosphere. Journal of Geophysical Research: Atmospheres, 124, 11070– 11085. https://doi.org/10.1029/2019JD030735
Zhang, W., Chen, G., Zhang, S., Gong, W., Chen, F., He, Z., et al. (2020). Statistical study of the midlatitude mesospheric vertical winds observed by the Wuhan and Beijing MST radars in China. Journal of Geophysical Research: Atmospheres, 125, 2020JD032776. https://doi.org/10.1029/2020JD032776

CC5: 'Reply on CC4', Gunter Stober, 25 Jun 2021
Dear Ryan,
to avoid misreading of my comment on the WGS84 geometry. We had a scientific presentation and agreed to share all the geometric aspects necessary for the retrievals. Please don’t consider it as an accusation. This was not intended. So my apologies if it is someone understood like that.
I suggest it is best to cite Clahsen (2018) and Stober et al., 2018.
Best regards,
Gunter

AC5: 'Reply on CC4', Ryan Volz, 25 Jun 2021
1. (Comparison to classical 2DVAR data assimilation) Sorry for the confusion, but we don't mean to say that the estimates are inherently done through 2D slices. The Gaussian processes are 4D, and the estimates at any point rely on all of the measurements in the manner described by the 4D covariance functions. In stating a 3 km vertical resolution, we mean to say that the covariance function we have used has a vertical length scale of 3 km, so that for any given estimate the measurements that contribute most strongly are within that vertical window. In this sense, one can think of the covariance function as having the effect of a type of soft binning.
2. (WGS84 geometry) We will be adding the Clahsen (2018) and Stober et al. (2018) citations in reference to the meteor processing in the revised manuscript, because the WGS84 aspect of the processing detailed therein is important enough to call out explicitly and not just transitively through Chau and Clahsen (2019). We appreciate the clarification.
3. (Vertical winds) The validity of the shortterm vertical wind estimates is still an open issue; what we claim is that our method provides uncertainties on the estimates based on the geometry and error propagation. We appreciate the references provided, but all of them deal with climatological vertical velocity estimates. In any case, we will leave the validation of such vertical velocity estimates for a future effort.

CC6: 'Reply on AC5', Gunter Stober, 28 Jun 2021
Dear Ryan,
thank you for the reply. The classification as 2DVAR refers here to the projection shown in equation 1. There is no explicit dependency in the vertical coordinate and as they mentioned in the first reply the retrieval doesn’t include vertical shears. For me, it appears that the presented retrieval could be used to obtain 3D winds in a 2D layer at a given time and altitude without knowledge of the layer below or above as well as the time step and before and after. That’s why I would classify it as 2DVAR. The covariance matrix is always 4D for such retrievals as we obtain solutions at the x,y,z spatial coordinates for different times. However, it is often more computationally efficient to solve blockdiagonal covariance matrices with just softwaredefined weighting. In so far, a 4D covariance does not present a unique token of the retrieval compared to retrievals for such applications.
Vertical winds are important. In particular, the magnitude of the winds is under debate. The references that were provided included climatological variability. In Figure 7. the manuscript presents winds with magnitudes of about +/ 510 m/s at scales of hundred kilometers for 15 minutes resolution using a 30 averaging windows (see covariance matrix). Such patches can be found in all panels. Considering the climatological knowledge and other observations that could be found in the references attached to the first comment, this discrepancy seems to be rather large and deserves some reflection in the discussion section. These winds seem to exceed by 200400% the magnitudes of other observations and model predictions.
It is also surprising that the paper references several times Vierinen et al. 2019. This paper basically contradicts the need for such more complicated retrievals. The basic claim in Vierinen et a., 2019 is that only radial correlations are needed to obtain Reynolds stresses without any inversion of the 3D winds on spatial grids. Vierinen et al., 2019 seem to be only applicable to the small subset of meteor detections, which are observed at the same spatial coordinates in x,y,z, and t in multiple links (at least 3). Such a discussion would strengthen the value of the presented manuscript and put more emphasis on why it is beneficial to apply more sophisticated mathematical approaches.
I am looking forward to the final published version.
Best regards,
Gunter

AC8: 'Reply on CC6', Ryan Volz, 16 Sep 2021
Thanks for the contributions from this comment thread. The discussion here, along with related comments by the other reviewers, will feed into changes in the revised manuscript.
In regards to Vierinen et al. (2019), we think that GPR and the correlation function technique described in that paper are complementary, giving different views to the same information. Analyzing just the secondorder statistics in that manner can give insight into how to improve the covariance model used for GPR. Extracting the full wind fields from GPR can give additional context to the higherorder measurements, like momentum transfer, produced by the correlation functions.

AC8: 'Reply on CC6', Ryan Volz, 16 Sep 2021

CC6: 'Reply on AC5', Gunter Stober, 28 Jun 2021

CC5: 'Reply on CC4', Gunter Stober, 25 Jun 2021

CC7: 'Reply on AC4', Aaron Smith, 29 Jun 2021
There seem to be a controversy regarding the prior use of WGS84 geometry. The use of WGS84 ellipsoid geometry is clearly described in the Stober et al. 2018 manuscript: https://doi.org/10.5194/amt1148912018. However, the cited M.S. thesis from M. Clahsen does not have an ISBN/DOI code and thus is not available for readers. I found the website of LeibnizInstitute of Atmospheric Physics at Rostock University providing their thesis works:
https://www.iapkborn.de/en/research/publications/thesismasterdiplomaphd
Unfortunately, no M.S. thesis from M. Clahsen is listed. There is a bachelor thesis from M. Clahsen (in German), dating back to 2015, which does not mention the WGS84 geometry. It would be beneficial for readers to share the cited M.S. thesis, either through the Institute website, or, if not possible, in the discussion forum here.

CC4: 'Reply on AC4', Gunter Stober, 25 Jun 2021

AC4: 'Reply on CC3', Ryan Volz, 24 Jun 2021

RC2: 'Comment on amt202140', Anonymous Referee #2, 18 Aug 2021
The present manuscript describes a new algorithm to estimate mesospheric and lower thermospheric (MLT) winds applying the Gaussian process regression (GPR) technique using meteor radar data collected during the SIMONe 2018 campaign in Germany. The algorithm is effectively a Bayesian approach, thus it requires prior information that is provided in the form of a mean and covariance matrices for the prior wind components. The covariance matrix is modeled as a Matern5/2 covariance function that has amplitude and scaling parameters that can be estimated as part of this technique or that can be provided directly based on some previous analysis of the observations. The document is well written, the description of the algorithm is clear, and the analysis of the results is profuse. Given the novelty of this approach and the need of measuring MLT winds with higher spatial and temporal resolutions to study the dynamics of this region, I would like to recomend the present manuscript to be published after some minor revision. Below, I present a few questions or comments that should be addressed by the authors.
* Based on the description of the GPR algorithm to estimate winds at any particular point in space and time, it is necessary to provide the mean and covariance matrix of the a priori wind distribution. While there is a nice description about how the covariance matrix can be modeled as a Matern5/2 covariance function, there is no much discussion about how to determine the prior mean wind. In the manuscript, it is mentioned that it is obtained from applying a tensor product cubic spline to the dataset. However, it is not clear, at least to me, if the authors are calculating this mean from wind estimates that were obtained applying a different method, for instance with the zeroorder method. Please clarify this point. In addition, the role of the prior covariance is analyzed in detailed, but the role of the mean is not. I can imagine that the results have a strong dependence on the a priori mean that is provided. The errors will probably increase depending on the accuracy of the mean. I would also suggest to discuss about the role of the a priori mean wind field in the manuscript.
* The position of the detected meteors are also the result of a fitting procedure, thus there are uncertainties associated to the space and time location of a meteor. How these uncertainties are taking into account in the GPR algorithm? Do they play any role in the accuracy of the estimated winds? I understand that a filter criteria is applied to the data to consider only high quality detections, what is the criteria that is used? Do the high quality detections have negligible uncertainties for the meteor locations? I would also suggest to discuss this question in the manuscript.
* The detection of meteors in time can be modeled as a Poisson process, in the sense that given a detection the probability of detecting the next meteor increases as function of time. Given this, the time location of the Doppler samples is also a Poisson process. However, for the GPR approach, we are assuming that the Doppler samples can be modeled as a Gaussian process in space and also in time. What is the impact of this difference in the estimation of the MLT winds?
* While the role of the covariance amplitudes for the apriori wind distribution are discussed in detailed, there is no much discussion about the role of the space and time scales. In principle, these parameters are also obtained from minimizing the negative loglikelihood, however, I can imagine that their values will strongly depend on the distribution of data considered. For instance given some particular data set where data samples are more concentrated around 90 km but more sparse around 80 or 100 km altitude, it is reasonable to expect that the space scale values will also vary as function of height, they will be probably shorter round 90 km but longer at higher or lower altitudes. Is this something that is considered in the algorithm? How the wind estimations will be affected by the precision of the space scales used in the algorithm?
* In section 5.2, in the case of the horizontal wind, it is not clear whether the “mean bias error” or the “mean absolute error” were calculated. The labels in Figure 3 indicates “u+v”, does this mean that the values of the zonal and meridional winds were added before calculating the error? Please, clarify this issue, the authors may consider to include the actual formulas used to estimate the errors. Also, in the same figure, the titles of the plots on the right side indicate “variance of horizontal wind”, however, I think the authors are referring to the variance of the horizontal wind error. Please, fix the labels of the figure to clarify their meaning. Similar comment with respect to Figure 4, have the authors plotted the variance of the vertical wind or the variance of the error?
* In section 6, in the implementation of the algorithm, it is mentioned that the data is divided in time intervals of 90 minutes with 30 minutes overlapping, would not this have an impact on the smoothness of the winds that are estimated? Would not it have sense to apply a similar criteria in space (altitude, latitude, and longitude) based on the actual distribution of meteors? The accuracy of the wind estimates in the areas further from the center may improve if covariance parameters are computed particularly for these regions. This comment is related to the role of the space and time scales presented above.
* In Figure 9, Gradiente winds and GPR estimates are compared. In particular, it is mentioned that the GPR estimates show some mesoscale structure. Are the authors implying that the GPR method was capable of estimating these mesoscale structure while the other method could not do it? I would recommend to clarify this comparison. The gradient winds were computed considering wider time intervals, and probably that has an impact in the smoothness of the estimated winds, thus the comparison would not be fair.
* In Figure 2 and Figure 7, the vertical winds are depicted as the color of the horizontal wind lines. In fact in the labels, it is mentioned “vertical wind speed”, however, “speed” by definition is a positive value, and the colors go from positive to negative, so I think it is more appropriate to change the label to “vertical wind”. In addition, are the values of the estimated vertical winds within the expectation? How do they compare with their corresponding variances? Are the error bounds small enough to have good vertical wind estimates? It would be important to add a discussion in the manuscript about the accuracy of the vertical winds estimated with the GPR method given that previous methods have just assumed that the vertical winds are zero.
* Finally, in the conclusion section, the authors indicate that the GPR method can resolved winds at the finest spatial and temporal scales allowable by the instrument. However, what are these finest scales allowable by the instrument? The geometry and the spatetime distribution of meteors will definitely set limits to the features that can be resolved both in space and time with this method, however, aren’t other methods within their assumptions also capable of resolving fine structures? In fact, in section 6, it is mentioned that the authors do not have a ground truth to validate the horizontal scales that are resolved by the GPR method. Given this it cannot be claimed that the GPR method resolves the finest scales allowable. I would recommend to change this conclusion.

AC6: 'Reply on RC2', Ryan Volz, 16 Sep 2021
Thanks for your detailed comments!
1. (More discussion of the prior wind mean) We concentrated on specifying and analyzing the prior covariance function because it is by far the more important component in the GPR specification, but we agree that more can and should be said about the prior mean function. The biggest effect of specifying a good mean function is that it improves the covariance function specification, allowing the amplitude and length scale hyperparameters to be smaller in general. This leads to smaller error bars on the wind estimates, but the wind estimates themselves (perhaps surprisingly?) don't change too much. We will be adding more discussion to the revised manuscript about how we specified the mean and its implications to the final estimates.
2. (Uncertainty in the meteor coordinates) Our current GPR method does not incorporate uncertainty in the meteor coordinates (space and time). We agree that it would be great to include this, but the GPR framework does not naturally incorporate this and adding it would be a significant project for future work. We expect this would entail leaving the closedform solutions behind and numerically sampling from the distributions (e.g. MCMC). But for the analysis in this paper, we do try to limit the effect of coordinate uncertainty by throwing out lowquality detections. We will clarify what this entails in the revised manuscript. Overall the uncertainties for the high quality detections are small enough relative to the covariance length scales that the added estimation error is negligible.
3. (Poisson process for meteor occurrence in time) One must be careful to distinguish between the probability distribution of meteors occurrence/detection and the distributions used to model the winds. For the wind process, when and where the meteors occur is irrelevant; all we care about is that the meteors produce a set of measurements, and how those measurements are distributed does not factor into our assumptions about the wind processes. A practical effect of the Poisson distribution for meteor detections, however, does mean that our coordinate sampling of the winds is more grouped than it would be if the meteors were uniformly random. That just means that we'll get lower uncertainties for the wind estimates in those regions due to the abundance of samples. It may also be relevant to point out that we will be adding more discussion of the wind process distribution assumption in response to other reviewer comments.
4. (Further analysis of the covariance length scale hyperparameters) This comment matches our experience with the length scales: the fitted values are strongly driven by the density of meteor sampling within a particular dataset, so we might naturally want to use smaller values around 90 km and during the morning detection peak and larger values at low/high altitudes and during the evening detection valley. This is not currently considered in our GPR technique, since we are using constant values for the length scales that don't vary with location. We think that allowing the length scales to vary with location (e.g. altitude) would likely lead to better wind estimates, and we think that this would be a fruitful area for future work (and have already suggested this in the manuscript, if not quite as clearly).
5. (Confusion on mean error, figures 3 and 4) We intend to say "magnitude of the mean error of the horizontal wind vector" (averaged difference in the horizontal wind vector, including both u and v components) for the bias plot, and we intend to say "variance of the horizontal wind *error*" for the variance plot. Same for the vertical winds. We will be explicit and update the manuscript/figures to refer to the "magnitude of mean vector error" or "mean error" and "variance of the error" or "error variance".
6. (Computation in overlapping intervals) The overlapping estimation procedure is not a necessary component of GPR, but it is helpful for reducing the computational burden as long as care is taken. And by that we mean that we have only made estimates at times when at least 45 minutes of data both before and after are included, for a total time window of 90 minutes (with centered estimate). This window is wide enough, given the 15 minute time length scale for the covariance, to ensure that the estimates produced are only negligibly different from the result of if a wider time window (or the whole dataset) was used. We have verified this by comparing the 90minutewindow estimates to ones done with a 180 minute window. Thus, the smoothness of the wind estimates is not affected.
It is an astute observation that the overlapping estimation procedure could be used to apply different hyperparameters that are more tuned to different segments of the data. This is indeed the most straightforward way to apply that type of analysis for future work, even if it is not terribly elegant. This is also how we know that the density of meteor detections affects the hyperparameters: we have observed that the fitted length scales in particular change somewhat throughout the day (when fitted to these overlapping windows of data), seemingly in correlation to the density of meteor detections. Fortunately the estimates are not changed greatly by imposing constant conservative values throughout the day; it just means that we're not achieving quite the best resolution at times when the meteor density would support it, effectively smoothing over potentiallydetectable features.
7. (Comparison to gradient method and estimation of mesoscale structures) Perhaps it is not totally fair to make this comparison and the claim that GPR shows mesoscale structure where the gradient method does not. It is definitely possible to perform an analysis with the gradient method (or other existing methods) that focuses on time and length scales similar to the GPR method, and thereby likely identify the same mesoscale structures. Such information is in the data, and we don't mean to claim that GPR performs some magic that unlocks it that is inaccessible to other methods. The benefit of GPR is not necessarily that it allows one to see these mesoscale structures, but that it provides a suitable framework and procedure for identifying those scales within the data and making them clear without manual data analysis.
8. (Vertical winds) In response to this and other discussion of the vertical winds and the figures, we have decided to remove the vertical wind component from the Figures 2 and 7 to improve clarity. Nevertheless, we will be adding more discussion of the vertical winds to the manuscript to address the questions raised in this and other reviews. The basic conclusion is that the technique is agnostic to the prior assumptions the user wants to employ for the vertical winds, and it also provides the necessary uncertainty information on the wind estimates that will allow the user to assess the quality of the vertical wind estimates. Through the typical meteor observation geometries, there is much less information about winds in the vertical direction than the horizontal directions. The fitting procedure on the SIMONe dataset produced a prior variance for the vertical wind component of about 90 m^2/s^2 using a set mean of zero. This could be from actual instantaneous nonzero vertical wind values, but it could also be elevated due to errors in the Bragg vector direction and/or meteor location causing contamination from the horizontal winds. The values produced in the estimates conform to this prior distribution and the information added through the measurements, but the posterior error bars are still large enough that a zero or nearlyzero vertical wind is a plausible explanation, especially considering the possible role of horizontal contamination. Great care is still needed in this and any future analysis of vertical winds, but we think GPR will provide a useful new tool in performing that analysis.
9. (GPR resolution) The discussion of comment (7) is also relevant here. It is evident that we need to clarify the point we are trying to make with this discussion and conclusion. We will update the manuscript to better highlight the ease with which GPR enables analysis at finer scales.

AC6: 'Reply on RC2', Ryan Volz, 16 Sep 2021
Ryan Volz et al.
Ryan Volz et al.
Viewed
HTML  XML  Total  Supplement  BibTeX  EndNote  

620  251  40  911  42  6  5 
 HTML: 620
 PDF: 251
 XML: 40
 Total: 911
 Supplement: 42
 BibTeX: 6
 EndNote: 5
Viewed (geographical distribution)
Country  #  Views  % 

Total:  0 
HTML:  0 
PDF:  0 
XML:  0 
 1