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 "non-parametric") We describe GPR as non-parametric 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 finite-dimensional 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 non-parametric 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 Matern-5/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 non-gridded 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 inter-comparison 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)?

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 second-order 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 cross-validation) Generalized cross-validation 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 cross-validation 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.

We appreciate your concern. We believe there is a misunderstanding, since our Doppler-derived line-of-sight 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.

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 well-known. 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 well-characterized 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 References

Chau, J. L., Urco, J. M., Avsarkisov, V., Vierinen, J. P., Latteck, R., Hall, C. M., & Tsutsumi, M. (2020). Four-Dimensional Quantification of Kelvin-Helmholtz 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 high-to-medium frequency mesoscale gravity waves by lidar-measured vertical winds and temperatures in the MLT. Journal of Atmospheric and Solar-Terrestrial 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 Solar-Terrestrial Physics, 105–106, 253–261. https://doi.org/10.1016/j.jastp.2013.08.012

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 closed-form 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 low-quality 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 90-minute-window 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 potentially-detectable 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 non-zero 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 nearly-zero 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.