Assimilation of atmospheric infrasound data to constrain tropospheric and stratospheric winds

This data assimilation study exploits infrasound from explosions to probe an atmospheric wind component from the ground up to stratospheric altitudes. Planned explosions of old ammunition in Finland generate transient infrasound waves that travel through the atmosphere. These waves are partially reflected back towards the ground from stratospheric levels, and are detected at a receiver station located in northern Norway at 178 km almost due North from the explosion site. The difference between the true horizontal direction towards the source and the backazimuth direction (the horizontal direction of arrival) of the incoming infrasound wave-fronts, in combination with the pulse propagation time, are exploited to provide an estimate of the average cross-wind component in the penetrated atmosphere. We perform offline assimilation experiments with an ensemble Kalman filter and these observations, using the ERA5 ensemble reanalysis atmospheric product as background (prior) for the wind at different vertical levels. We demonstrate that information from both sources can be combined to obtain analysis (posterior) estimates of cross-winds at different vertical levels of the atmospheric slice between the explosion site and the recording station. The assimilation makes greatest impact at the 12-60 km levels, with some changes with respect to the prior of the order of 0.1-1.0 m/s, which is a magnitude larger than the typical standard deviation of the ERA5 background. The reduction of background variance in the higher levels often reached 2-5%. This is the first published study demonstrating techniques to implement assimilation of infrasound data into atmospheric models. It paves the way for further exploration in the use of infrasound observations - especially natural and continuous sources - to probe the middle atmospheric dynamics and to assimilate these data into atmospheric model products.


I. INTRODUCTION
Despite much recent attention to extra-tropical stratospheric dynamics and their connection to the troposphere, the amount of observational data in the stratosphere available to numerical weather prediction centres remains limited. A better representation of the stratospheric dynamics and the stratosphere-troposphere coupling in models has the potential to enhance tropospheric weather forecasts, in particular, on sub-seasonal timescales [Baldwin et al. 2003, Charlton and Polvani 2007, Haase et al. 2018, Karpechko et al. 2016, Kawatani et al. 2019, Kidston et al. 2015, Mitchell et al. 2013, Pedatella et al. 2018, Polavarapu et al. 2005, Taguchi 2018]. Moreover, the lid of several atmospheric model products has been raised into the mesosphere [Polavarapu et al. 2005] and it has been demonstrated that this can improve numerical weather and climate models , Kidston et al. 2015, Orsolini et al. 2011. But the full potential of high-top models can only be unlocked if middle atmospheric winds are better represented [Baker et al. 2014, Korhonen et al. 2019, Lee et al. 2019. Hence, it is timely to explore novel datasets and assimilation approaches that can constrain the upper stratospheric dynamics in atmospheric model products.
Infrasound waves are acoustic waves at frequencies below the human hearing limit (typically around 20 Hz). These waves can be generated by natural sources, such as volcanoes, earthquakes and ocean swell, but also by human sources, such as mining and explosions (see, e.g. Le Pichon et al. [2018]). These waves propagate through the atmosphere and can be recorded by ground-based stations. The wave frequencies of greatest interest for atmospheric characterisation are typically of the order of 1 Hz. The time and form of the received signals provide temperature and wind related information about the atmosphere the waves traverse. Infrasound waves may travel from sources on the surface of the Earth, reach a maximum altitude where they are partly or fully reflected or refracted, and then reach back to the surface to be detected by a receiver. Effectively, they probe a slab of the atmosphere in a tomographic fashion since the time it takes for these waves to complete their path is affected by the characteristics of the atmosphere they pass through: in particular, the wind velocity and temperature, but also attenuation-related properties like density and relative humidity. Hence, spatio-temporally integrated information carried by the propagating infrasound waves can be utilised to reconstruct or constrain atmospheric variables. Sound waves are already exploited in other tomographic and imaging problems. For instance, in underwater acoustics, temperature profiles [Dzieciuch et al. 2013] and seafloor bathymetries [Wölfl et al. 2019] are mapped using sound waves. Probabilistic infrasound propagation has been studied in Smets et al. [2015], where measured infrasound wavefront parameters for one year of infrasound explosions were compared to raytracing simulations using the ensemble atmospheric wind and temperature fields of the ECMWF ensemble data assimilation system of perturbed analyses [Buizza et al. 1999].
The current study follows directly from a recent paper by Blixt et al. [2019], which used the same dataset to demonstrate that atmospheric cross-winds can be estimated directly from infrasound data using propagation time and backazimuth deviation observations, and interpreted these results in the context of ERA-Interim reanalysis winds. There is a physical effect which is the basis of this work: when a steady cross-wind acts on a propagating acoustical plane wave, a bending of the wave-front is introduced. This creates a deviation in the apparent backazimuth direction of infrasound wave-fronts impinging on ground-based sensor array stations. We use this physical effect to assess the dynamical evolution of the stratosphere during several events, as sampled by the infrasound waves on their paths between Finland and a ground-based station in Northern Norway. The array signal processing algorithms exploit infrasound signals recorded on a set of 25 sensors distributed on the ground within a 3 km wide aperture (see Blixt et al. [2019], Figure 1).
Data assimilation (DA, e.g. Asch et al. [2016], Kalnay [2003]) is a discipline which aims to combine different imperfect and incomplete sources of information to produce a better estimate of a variable of interest. In particular, it takes into account the uncertainty of the information sources. The most ambitious approach obtains and updates descriptions of a system using probability distribution functions (pdf's) by application of Bayes' theorem. In practice, however, sample estimators like mean, covariance and mode of the distributions often suffice. In particular, the Kalman filter [Kalman 1960, Kalman andBucy 1961] and its ensemble implementation [Burgers et al. 1998, Evensen 1994, Tippett et al. 2003] assume Gaussian statistics in the sources of errors, as well as none or small deviations from linearity in the evolution and observation processes. The filter operates with the first two statistical moments of a distribution. An advantageous feature of the Kalman filter is that it can assimilate an integrated observation variable (in our case an average wind component resulting from vertical integration along the path of propagation) and translate this into increments at different vertical levels. This proved useful, for instance, in the assimilation of radiance satellite observations [Lei et al. 2018]. A discussion on the prospects of assimilating atmospheric infrasound data into numerical weather prediction models can be found in Assink et al. [2019].
There are two main objectives of this study. The first is to develop a framework which allows for assimilation of tropospheric and stratospheric wind information based on atmospheric infrasound data. The second is to provide a first demonstration and proof-of-concept with an offline (i.e. no cycling involved) infrasound DA experiment using the developed framework, exploiting a dataset which is already well-characterised in previous works.
We generate an estimate of the averaged cross-wind component along the relevant track from the explosion site in Finland to the station in Northern Norway, as well as an associated measure of uncertainty. We apply the deterministic ensemble Kalman Filter (DEnKF) as described in Sakov and Oke [2008]. We select this approach because it allows for model-space localisation, as opposed to observation-space localisation which is not feasible for integrated quantities [Lei et al. 2018]. Some specifics of this method are outlined in Appendix A.
The remainder of this paper is organised as follows: Section II explains the system setup, detailing the way observations are related to the state variables of the system under different degrees of simplification from the most general problem to the case considered in the current work. In Section III, we perform synthetic-data experiments under ideal conditions with an infinite ensemble size, and with different vertical weights in the observation operator. These experiments verify the offline DA process in a controlled setting. Section IV presents the real-data assimilation experiments using infrasound from 18 years of explosions. In Section V we conclude the study, discuss its limitations, and provide ideas and suggestions for future work.

II. CROSS-WIND EFFECTS ON THE PROPAGATION AND ARRIVAL OF INFRASOUND WAVEFRONTS
Let us explore the effect of a cross-wind on the propagation of infrasound waves. Recall the basic principle: a background wind field affects the propagation of infrasound waves; specifically, a cross-wind can bend the wave front. Infrasound waves, however, do not modify the background wind field.

A. Propagation within a plane
First, we discuss horizontal wave propagation only. We illustrate the situation in the left panel of Figure 1. Consider a plane with two horizontal directions denoted r a and r c ; the sub-indices a and c denote along-track and cross respectively, and refer to the wind direction with respect to the propagation of the infrasound wave. The infrasound source is the red star in the bottom, labelled S, while the receiver is the red star in the top, labelled R. The straight red line connecting the two points has length d a . The backazimuth θ is the angle of this line measured with respect to the North. Now consider a constant (for now) wind W c blowing perpendicular to the line d a . The effect of W c is to create a change in the apparent backazimuth direction of infrasound wave fronts when they arrive to R [Diamond 1964]. The received waves seem to come from an apparent source marked by the blue star S , a distance d away from R, and with a modified backazimuth angle θ . The distance between the real and apparent sources (purple line) is d c . The change in angle is denoted as: (1) For positive cross-winds (as in the setup of the figure) the change in angle is negative. We can relate different elements of this system using the following considerations. The wave is emitted from S and received at a time T after the explosion. For infrasound waves propagating within atmospheric waveguides, celerity υ is defined as the ratio between the straight distance between source and receiver divided by travel time, i.e: The lines d a and d c are the two legs of a right-angled triangle. We can solve for d a from (2). For d c we simply have (since the cross-wind is constant): Some trigonometry yields tan (|∆Θ|) = d c /d a . Explicitly, this is: The negative sign comes from the direction ∆θ that is defined in (1). The middle panel of Figure 1 illustrates (4) for different values of cross wind (horizontal axis) and celerity (lines). Note that as long as W c υ the function is close to linear. This is verified by the McLaurin expansion of the arctangent function: So far, we have considered a constant cross-wind W c . In the general case this speed can be a function of the position r a and time t, i.e. w c (r a , t). Then (3) becomes an integral: This computation is not easy in general. The position r a of the wave front depends on the infrasound speed of propagation along d a , which is the sum of the sound speed (a function mainly of temperature) and the actual background wind along the direction of the propagation w a (r a , t). Blixt et al. [2019] define an average cross-wind as: With this definition we can still use (4), with W c being the average cross-wind velocity along d a . It is actually more useful to convert this time integral into a spatial sum. We illustrate this process with the aid of the diagram in the right of figure 1. We divide the line d a into N segments. Each segment has length d a n , and it is clear that: Similarly, the total travel time T is the sum of the time spent in each segment: Consider also a constant background sound speed C. Also, consider that the two wind components are constant per segment of d a . In the n th segment d a n we have: {W a n , W c n }. The time for the infrasound wave to travel a segment d a n is: Following (3), the cross displacement in the n th segment is the product: d c n = W c n T n . The total cross displacement is then the sum: Notice that by dividing the last expression by the total travel time T , we can define a weighted average cross wind speed as: with the weights: Then, we can still use (4) to relate this weighted average to the change in backazimuth angle. Most importantly, we can estimate the average cross wind as a spatially weighted linear combination of cross winds. The weights derived in (12) follow from several simplifications. There are wave-tracing techniques that can model the trajectory of the propagating wave (see, e.g. Hedlin and Walker [2013]) which can be used to determine adequate weights.

B. 3D Propagation
Having explained the basics, we now move to full 3D wave propagation, i.e. when the trajectory of the infrasound wave has a vertical component. This is depicted in Figure 2. The top panel shows an atmospheric volume discretised to a model grid. Both the source (S) and receiver (R) are in the surface. In this case, the line d a is a segment of the great circle between S and R, and it is not necessarily aligned with the grid. A simple example path of an infrasound wave is shown in yellow. The wave travels both in the r a and z directions. The wave travels in the vertical to a given maximum altitude from where it returns down to ground (e.g. due to partial reflection as explained in Blixt et al. [2019]) and it is then detected at the receiver. As in the 2D case, the wave travels through a cross-wind field which leads to a change on backazimuth angle ∆θ towards an apparent source S . This cross-wind now also depends on altitude: w c (r a , z, t).
As before, the travel time T is known, as well as the horizontal distance d a , so we can still define celerity as in (2). The expression for d c is the same as (6), but in this case the w c also depends on altitude: Turning this time integral into a spatial sum is slightly more complicated. The process is illustrated in the middle panel of figure 2. First, the situation is reduced to a 2D problem by creating a channel centred in the line d a . The winds from the native grid are interpolated to provide the along and cross values in this new setup. As before, we divide the distance d a into N different intervals, and this time the distance from z = 0 to z = z max (which has to be determined) is divided into N z intervals. This creates a 2-dimensional DA grid where we consider the {n, n z } th box to have a cross-wind W c n,nz , which is obtained from the native grid via interpolation. It is important to notice that the wave does not go through all the boxes, but only a set of them, referred as Ω below. The displacement d c is affected only by the cross-wind in the boxes of this valid set.
The expression for the average cross-wind, i.e. the equivalent to (12), becomes a double sum: The weight is zero for any box outside the set Ω. For the boxes in the set Ω the weights are more complicated than in (13), since the (diagonal) length travelled by the wave in different boxes may be different, and the effective propagation speed includes both the along-track wind and vertical velocities. Therefore, one may rely on ray-tracing techniques to derive these weights.
Since this is our first study, and the distance between source and receiver is relatively short (178 km), we further simplify the problem as illustrated in the bottom panel of figure 2. We do this by considering only N = 1 interval along the propagation of the wave. For the rest of the work we consider the average cross-wind speed as a weighted sum over N z vertical levels: In DA terminology, our state variable is the vector of vertical (horizontally averaged) cross-winds w c ∈ R Nz : It is also useful to group the vertical weights in a vector as: In this case the sum (16) is simply a vector product: and (4), in DA terms the observation equation, can simply be written as: This is mapping that reduces dimensionality R Nz → R.

III. SYNTHETIC-DATA ASSIMILATION EXPERIMENTS
This section describes basic synthetic DA experiments, before moving to the case of the assimilation of recorded infrasound data. Consider N z = 4 vertical levels in a propagation volume. Let the cross-wind w c ∈ R 4 be a Gaussian random variable with zero mean µ b = 0 and covariance B ∈ R 4×4 . We apply the DEnKF with sample size of N e = 10 4 elements. This sample size is practically free of sampling noise. This allows for the computation of accurate estimates of the associated pdf's, and for the background ensemble mean and covariance to be virtually identical to the real ones, i.e.x → µ b and P b → B.
The ensemble background covariance P b is crucial since it spreads information from observed to unobserved variables, see Appendix A for details. In Section IV, where we perform offline DA with real measurements, we estimate the background covariance from an ensemble reanalysis model product. In the current synthetic example, however, we prescribe a background error variance which is constant at all vertical levels, such that B can simply be written as a product of a common variance (scalar) and a correlation matrix: We set the variance to σ b 2 = (10 m/s) 2 . We prescribe two correlation matrices C ∈ R Nz×Nz with only positive correlations. The ij th elements of these matrices are: In the first case δ i,j is the Kronecker delta, so C becomes the identity matrix. The second case renders a Toeplitz matrix with a main diagonal of ones and an exponential decay for the off-diagonal elements. Both matrices are plotted in Figure 3 for visualisation. We also prescribe the altitude-dependent weights α ∈ R 4 applied in (17). We consider three cases: For α low , the effective cross-wind simply becomes the cross-wind at the lowermost level, while for α all , the effective cross-wind becomes the averaged cross-wind over all four altitude levels. For α top , the effective cross-wind is the average of the two cross-winds at the highest levels. This case is less realistic, but included for comparison. The celerity is set to the fixed value υ = 300 m/s, and we assimilate an observation with a given value and the prescribed uncertainty: Solving from (20) yields w c ≈ 60± 6.2 m/s. We find this approximate corresponding error using the linear approximation: Figure 4 shows the results of the assimilation experiments considering the two matrices P b given above. This figure has 3 columns, one for each set of weights α. We plot several pdf's in each panel. To ease visualisation, the pdf's are scaled, and hence the vertical axes have no units. The background pdf estimated from the model ensemble is shown with a grey dotted line for the four vertical levels and operators. We also plot the analysis pdf's for the two covariances. When P b is diagonal, the DA process can only update the levels with non-zero values in α. The analysis pdf's corresponding to this case are shown by black dashed lines. In the left column, only the lowest level is updated, while in the centre column the 4 levels are updated. In the right column, only the top two levels are updated. All observed levels are updated similarly as we apply a non-zero operator with equal values.
A non-diagonal covariance matrix P b yields a different result because non-zero off-diagonal values communicate information from observed to unobserved levels. The blue dotted lines show the analysis pdf's for this case. The magnitude of the update decreases with distance between the observed layers and non-observed layers, as expected from the exponential off-diagonal decay in P b .

IV. OFFLINE ASSIMILATION EXPERIMENTS USING OBSERVED INFRASOUND FROM EXPLOSIONS
We finally proceed to perform offline DA based on real infrasound recordings. The offline character implies that the assimilation at a given observed time is independent from all other times.

A. Observations and background
Our observations come from a dataset recording explosions at the Hukkakero site in northern Finland [Blixt et al. 2019. These explosions series are conducted during August and September, with individual explosions typically separated by about 24 hours. The dataset considered in the current study covers the years 2001-2018. The infrasound waves produced by these explosions are detected at the ground-based ARCES array station in Norway, which is located 178 km due north from the explosion site. It takes the wave around 10 minutes (in average) to propagate from the source to the station.
Since we know the exact explosion and detection times, as well as the exact source and receiver locations, we can compute the celerity υ value with high accuracy. In fact, we will consider it to be error-free. The backazimuth deviation angle ∆θ for each explosion is obtained from observations. For these observations we consider an unbiased error following a normal distribution with a standard deviation of 1/20 of a degree. See Blixt et al. [2019] or, e.g., Szuberla and Olson [2004] for details on the estimation of observational error in this case. Figure 5 displays the backazimuth deviation ∆θ (top panel) and the celerity υ (bottom panel) for each explosion. The years are separated by black vertical lines. To facilitate visualization we do not display the exact time of each explosion. We discard data points where the magnitude of the backazimuth deviation is |∆θ| ≥ 0.75 rad (not shown in the figure), retaining a total of N = 370 valid events. Table I lists the number of events used and discarded for each summer.
We extract the background cross-winds from the ERA5 reanalysis product [Hersbach et al. 2019], which has 10 ensemble members. We interpolate the horizontal winds from the native grid to the along-track and cross directions to the great circle connecting Hukkakero and ARCES. This is done for all the 137 ERA5 vertical levels. The time resolution of ERA5 ensemble product is 3 hours, so we linearly interpolate the wind values to the origin time of the explosion. The propagation time from source to receiver, which is around 10 minutes, is disregarded when extracting the ERA5 winds. This simplification would not be valid for longer propagation times. Figure 6 shows statistics for the background cross-wind velocities for the 137 vertical levels (vertical axis) at the time of each explosion over the 18 years (horizontal axis). The vertical lines show the change of year and again the exact times are not shown in the axis. Note that the vertical levels do not have uniform resolution. The top panel displays the sample mean over the 10 ensemble members. We scale the colours to cover W c nz ∈ [−25, 25] m/s. In general, the mean cross-wind speed is characterised by a strong positive jet in the lower levels (around z = 10 km), and a strong negative cross-wind in the upper levels (around z = 60 km). The cross-wind shows, however, a significant variation in time.
The bottom panel of Figure 6 shows the cross-wind sample standard deviation over the 10 ensemble members. Lower levels have smaller standard deviations than higher levels. For instance, the region above z = 50 km has standard deviations of up to 2 m/s or larger, whereas the standard deviation in levels below 30 km are rarely higher than 0.5 m/s. This is expected since the reanalysis data contains information from atmospheric wind observations from these altitudes. The number of observations generally reduces with height [Duruisseau et al. 2017]. This plot suggests the observational impact of the infrasound measurements to be higher in the levels above around z = 30 km. However, the other factor for this impact involves the coefficients for different vertical levels, which is something we discuss in the next subsection. Figure 6 displays a time-varying black line at around z = 40 km. This represents the estimated maximum altitude the infrasound penetrates before being reflected towards ground. Any altitudes above those lines cannot be updated directly from the observations. Therefore, updates above this line are due to vertical covariances in the DA process. The return altitude of the infrasound is estimated by matching the travel time of a modelled infrasound ray through the model atmosphere with the observed infrasound travel time, as explained in Blixt et al. [2019].

B. Vertical weights
In the synthetic experiments we prescribed coefficients to compute the weighted cross-wind average. In the current section, we estimate these weights from ray-tracing through wind and temperatures [Blixt et al. 2019] extracted from the ERA-Interim reanalysis atmospheric product [Dee et al. 2011]. This is shown in Figure 7 for 14 events in 2016. The lines are coloured according to the corresponding celerity υ as indicated in the label box on the right of the figure.
This figure shows un-normalised vertical weightsα nz for each explosion. Notice that none of the explosion-generated infrasound waves penetrate higher than 50 km altitude, with the majority only reaching around 40 km. It is clear that the waves spend a significant part of the propagation time within the lowermost 10 km levels and within 30 and 40 km. The celerity υ ranges between υ = 274.4 m/s and υ = 292.9 m/s for these events.
This process is applied to all 370 explosions, yielding vertical coefficients and maximum vertical penetration values for all the events for 18 years. These profiles are plotted in Figure 8. The horizontal axis corresponds to time, the vertical axis to altitude, and the colours correspond to the un-normalised coefficients.

C. The data assimilation
To perform DA, we first need to define the vertical levels to use in the process. If we estimate the cross-wind at each reanalysis level, the size of the state variable becomes N x = 137. This problem is quite challenging, especially since we only have N y = 1 observation containing integrated information. An extra complication comes from our relatively small ensemble size (N e = 10). A way to simplify the problem is to create fewer DA vertical levels by applying vertical averaging. After trying several averaging kernel heights, we decided to use N zDA = 6 equidistant DA levels with a height of ∆Z DA = 12 km, covering the altitudes z = 0 km to z = 72 km. In a given DA level l with N zl non-equidistant reanalysis levels inside it, the vertically averaged cross-wind is: We use this weighted approach also to obtain a weight α l for each DA level. We also ensure the sum of the weights to be normalised one. Starting from the un-normalised weightsα at native levels coming from the ray tracing, we compute: The normalised weights computed using (26) are plotted in Figure 9 for all the events (horizontal axis) and each DA vertical level. There is temporal variability in the weights, especially for the lowermost four levels. Note that the uppermost level (60-72 km) always has zero weights, and in the next level (48-60 km) infrasound waves penetrate for only few events per year. These upper levels can only be affected by observations through the sample covariance between different levels.

D. The quality of the background covariance
An accurate representation of the background covariance matrix is vital to the DA process. Recalling that we have a limited-size ensemble, N e = 10, a low-quality estimator can be harmful to the analysis values. Localisation can handle poor-quality long-distance covariances when working with small ensemble sizes. In fact, even after reducing the problem to N zDA = 6 DA vertical levels, we are still left with noisy background error covariance matrices P b ∈ R 6 . This is illustrated in the top row of Figure 10, which shows the raw correlation matrices from the ensemble at 4 different times (columns). Green colours are positive correlations and purple colours are negative correlations. These matrices are different because they contain flow-dependent information. We display correlations and not covariances because the variances for different altitudes have different orders of magnitudes. Intuitively, the correlations should decrease with increasing vertical distance. The second row shows the covariances after being localised, which means they have been multiplied (using Schur product denoted •) times a matrix Φ of coefficients that decay with distance [Hamill et al. 2001]: We apply a Gaspari-Cohn localisation function -which is a compact-support approximation to a Gaussian [Gaspari and Cohn 1999] -with a half-width of 15 km. We choose this to be larger than the height of each DA vertical level (12km). Another technique to improve sample covariances is inflation [Anderson and Anderson 1999]. In its simplest implementation the ensemble of background perturbations is scaled: Inflating the background covariance increases the Kalman gain, which makes the observations have a larger impact in the analysis field. This makes intuitive sense: if the uncertainty in the background is considerably larger than the uncertainty in the observations, the assimilation should tend to ignore the background. There exist more advanced inflation implementations, e.g. where time and space-dependent coefficients are applied [Miyoshi 2011, but in the current experiment these are fixed. Moreover, it is common to choose a ρ value which minimises a set of accuracy measures -e.g. the root mean squared error of the analysis mean -with respect to independent observations. In the current work, we do not have reliable estimates of such verification values, so instead we study the impact of different inflation values. A clear reason for applying inflation is that the ensemble background covariance is often underestimated. This is inherent to small ensemble sizes, see Amezcua and van Leeuwen [2018], Sacher and Bartello [2008], van Leeuwen [1999] for detailed explanations of direct and indirect effects. There are more tangible mechanisms for the misrepresentation of the background covariance. This includes differences in the resolution of model and observations, and the imperfect representation in the forecast and observational process.
In our experiment setup we recognise there are sources of imperfection. These include (a) we temporally interpolate from the reanalysis times to the time of the observation, (b) we consider instantaneous velocities while the infrasound wave propagates for around 10 minutes, (c) there may be erroneous assumptions behind the calculation of the α weights inside the ray-tracing technique. We performed the experiments with several inflation values α, and below we discuss the results obtained using two of these values.

E. Results
Here, we display the results for the following DA settings: N x = N zDA = 6 state variables per observational instant, N y = 1 observations, N e = 10 ensemble members, vertical localisation with a half-width of 15 km, vertical weights coming from the ray-tracing assumed to be perfect, and two different inflation factors ρ = 0, ρ = 1. The second inflation value means the standard deviation of the background is doubled compared to the data in the non-inflated assimilation. Figure 11 shows the weighted cross-wind solved from (4) for observations (black line) and computed from (19) for the background (blue line), as well as the resulting analysis (red and green lines, depending on the inflation). To facilitate visualisation, we only display the years 2001 and 2002.
The background and observation cross-winds are similar for some of the events, but for most events the DA produces changes. In fact, for some events the difference is up to 1 or 2 m/s. In the absence of inflation, the background and analysis values are quite close. The use of inflation, however, increases the differences between analysis and background as expected.
The impact of the observations is in general low, especially in the absence of inflation. Several factors can explain this: First, the variance of the background ensemble is small, which is expected since this is a reanalysis product already containing information. The observation impact might be greater if instead using an ensemble forecast as background. Second, as already mentioned, the ensemble size is small with only N e = 10 members. A larger ensemble would allow to select different state variables, for instance a larger number of DA vertical levels. Less vertical averaging of the original variables would give a prior with larger variance, hence allowing for larger observational impact. It is important to point out that in an online setting, the background would come from an ensemble forecast and the infrasound observations would not be the only data assimilated. Another aspect is that the stratospheric winds are in general significantly weaker and less variable in August and September than in winter. It will be interesting, in a future work, to perform these experiments for wintertime explosions and to assess the observational impact.
How do changes in the vertically averaged cross-wind translate to the different DA vertical levels? These results are shown in Figure 12, which has two panels corresponding to two selected vertical levels: 0 − 12 km (bottom) and 48 − 60 km (top) for the 2001 and 2002 events. The blue line shows the background mean, with the cyan lines to each side indicating one standard deviation. The red line denotes the analysis mean, with the magenta lines to each side indicating one standard deviation. This analysis was produced using inflation. There are some changes in the values of the cross-wind in the lower level, however these tend to be small. The difference between background and analysis is more noticeable at higher DA levels, which are not even updated directly (recall most explosions do not penetrate these altitudes) but based on the inter-level covariances. In the no-inflation case, there are still changes, but they are less distinguishable in the plot.
To evaluate the results for all DA events without inflation, we compute the following derived quantities at each time t: The quantity (29a) is called analysis increment, which is the difference between the analysis mean and the background mean. The second quantity (29b) is a variance ratio, which is the analysis variance divided by the background variance. Both the analysis increment and the variance ratio are computed for each DA level and for each observation time.
These quantities are displayed in Figure 13 for all events as a function of time (horizontal axis) and altitude (vertical axis). The DA process does not only result in a modified mean, but it also reduces the uncertainty of the estimate. In mathematical terms, the trace of the analysis covariance matrix is smaller than the norm of the background covariance (see e.g. Asch et al. [2016]). The top panel displays the innovations with typical magnitudes between −0.2 m/s and 0.2 m/s. Pink colours represent positive increments and green colours represent negative increments. Remember that these increments are the changes that the infrasound observations produce to the forecast. The bottom panel displays the resulting variance ratios of standard deviations. The figure plot confirms that these values, as expected, are always smaller than 1. The darkest colours correspond to the greatest reduction in uncertainty. The experiment results in a ratio which descends to around 0.9. However, we keep in mind that the reanalysis data ensembles already contained small statistical uncertainty. Figure 14 summarises the increment d ab (left) and the variance ratios r ab (right) obtained, for each vertical level. Box plots provide a non-parametric summary (with outliers omitted to avoid cluttering the figure). These box plots are complemented by the mean as shown by blue dots. There are non-zero innovation results at all vertical levels, with the largest typically within the level 24 − 36 km, and the smallest typically within 48 − 60 km altitude. In the three upper-most altitude layers, at least 75% of the increments are negative. Note that in at least three levels, the mean and the median differ significantly, indicating asymmetry in the distribution of the innovations.
The right panel of Figure 14 shows the resulting variance ratio, which is 0 and 1. Since this has a non-symmetric distribution, the mean and median does not coincide. Notice that the reduction of the variance is largest in the four upper levels, i.e. 24 − 72 km. This is expected because these levels have greatest background uncertainty. Since the waves penetrate only to around 40 km, updating the DA levels at these altitudes is done both directly and via covariances. The lowermost two levels exhibit a limited covariance reduction. Although the coefficients for the lowest DA level are significant, the background winds are already well constrained there, hence only allowing the assimilated infrasound-based data to impact the analysis to a minor extent.

V. SUMMARY AND FUTURE WORK
This is the first study to explore assimilation of atmospheric infrasound data into atmospheric models in order to constrain atmospheric winds. The backazimuth deviation of infrasound waves carries integrated information related to the cross-winds acting on the wave along its atmospheric propagation path. We show that assimilation of this information using an ensemble Kalman Filter is able to provide corrections to the wind in stratospheric and tropospheric altitudes.
We performed DA experiments for 370 explosion events throughout 18 years (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018). We know the accurate time and location of the explosions and arrivals of the infrasound waves. This allows us to accurately calculate the propagation time and the celerity υ. It also permits performing complementary ray-tracing to determine the vertical sensitivity at different vertical levels, which is needed in the observation operator. To reduce the dimensionality of the problem, we consider average values corresponding to N z = 6 DA levels, each 12 km thick. This is opposed to the original N z = 137 levels of the reanalysis. Here, there might be room for improvements and subsequent works can explore in detail the effect of selecting different numbers of DA vertical levels.
The results of the DA experiments yielded non-zero analysis increments (defined as analysis mean minus background mean) for most times, with the largest values in the 24 − 36 km layer. More than 75% of the increments calculated above 36 km are negative, suggesting a bias in the background values. As required by construction, the variance in the cross-wind values at all levels has been reduced for all data points assimilated, while for the upper-most levels the reduction reaches up to 2 − 5%. This implies a reduction of the uncertainty in the estimation.
It would be desirable to apply this framework to existing datasets for explosions performed during the winter season when the stratosphere is more dynamic than in August and September. This may present a challenge, though, since larger magnitudes of cross-wind can reduce the linearity of arc-tangent in (4). This may present a challenge to the DEnKF. We can instead try with techniques that handle departures from linearity better. For instance, the iterative ensemble Kalman smoother (e.g. [Evensen et al. 2019] is a useful candidate to solve this problem.
For future work, we suggest exploiting signals from natural continuous sources like microbaroms. These are atmospheric infrasound waves produced by ocean surface hot-spots where counter-propagating surface waves are prevalent [De Carlo et al. 2020, den Ouden et al. 2020, Donn and Rind 1971, Le Pichon et al. 2006, Posmentier 1967.
In this work we rely on many simplifications. In the future we aim to solve a setup akin to the middle panel of Figure 2. Namely, this would consider the cross-wind variation both on the vertical and the along-track direction of the infrasound wave. This becomes especially important when the distance between source and receiver increases. An example is the detection in Norway of infrasound from microbaroms generated near Iceland; in this case the separation is about 2000 km and considering a single horizontal slab may be detrimental to the usefulness of the estimation. In this case we may also not be able to consider the winds constant in time for each position along the trajectory.
We have an important advantage when working with the explosions data set: the times and locations of both the emission and detection of the infrasound waves are known accurately. This, in turn, allows us to consider the celerity υ perfectly known, which we have done in this work. In the case of microbaroms, for instance, the time and location of the detection may be accurately known, but the location and time of the emission may prove much more elusive. In these cases, there may be large uncertainty on the values of celerity. This, added to the uncertainty in the propagation medium, makes it necessary to consider celerity as another random function. Several previous works establish the pdf for celerity in infrasound propagating under stratospheric waveguide conditions. For example, Blom et al. [2015] used simulations to establish the expected celerity to be between 250 and 350 m/s for propagation distances at around 200 km. Similarly, Morton and Arrowsmith [2014] analysed both simulations and measurements to find a celerity distribution at 275 km distance with values between around 280 and 310 m/s. A data-based study presented in Nippress et al. [2014] estimates the celerity distribution at 200 km distance to span the 270 to 300 m/s range.
Regarding the dynamics of the wave propagation, we recognise that the framework applied here requires an auxiliary ray-tracing method to determine the sensitivity to the wind at different vertical levels (first native and then DA levels) in the weighted sum giving the average cross-wind impacting the observation. Follow-up studies could include the development of approaches to instead estimate these sensitivity weights as part of the assimilation process. Then the implementation of an expression akin to (13) might be required. In turn, this would require the state variable to include the along-track wind and the temperature (which the sound speed is a function of) in each of the grid points traversed by the wave. However, this would also give an opportunity to estimate along-track winds.
An important detail to mention is that a DA process requires a verification step to assess the quality of the analysis field obtained. For identical-twin experiments, one produces a synthetic truth from which the simulated observations were extracted. Then the analysis can be assessed with respect to this reference truth. In operational DA the true state of the system is unknown, so verification becomes more elusive. One option is to have independent observations or independent reanalysis data that can verify the analysis. In the current study, we do not have independent observations for validation. This in turn restricted us from tuning the values of localisation radius and inflation parameter. Although this is outside the scope of the current study, prospective future studies might have access to independent measurements to allow for tuning and verification. Here, the ADM-Aeolus satellite mission will likely be a reliable benchmark for winds up to 30 km altitude [Tan et al. 2008]. Likewise, future validation may be possible using data from portable lidars; for example, the CORAL system [Kaifler et al. 2015[Kaifler et al. , 2017 might be upgraded to provide direct wind measurements.
Finally, the DA experiments of this study were made offline. In order to perform online assimilation experiments, it would be necessary to implement the methodology in a dynamic forecasting system. In an operational or quasioperational setting, infrasound measurements can be added to the rest of the available observations at the moment of assimilation. Although implementation in an operational assimilation system still requires substantial further work, the methodology described in the present study provides a starting point for such developments. This an objective of a next step following up the European ARISE and ARISE2 projects [Blanc et al. , 2019.
Given that single-station infrasound measurements provide atmospheric wind measurements within a sparsely observed altitude range for a given geographical region, an extended or even global multi-station wind sampling might be feasible using, e.g. infrasound station data recorded by the International Monitoring System network [Dahlman et al. 2009, Marty 2019. Hence, there are several opportunities yet to explore in further work related to atmospheric probing and data assimilation using infrasound datasets. A long-term objective is to enhance or constrain the representation of stratospheric winds in global models, thereby contributing to enhanced surface weather predictions on subseasonal-to-seasonal timescales [Domeisen et al. 2020a,b].  (S) and receiver (R) and the apparent source (S') and the receiver. The middle panel shows the relationship between change in backazimuth angle and the cross-wind for different celerity values. The right panel shows the same situation as the left panel, but after spatial discretisation. Figure 2. Depiction of an infrasound wave (yellow) travelling through an atmospheric volume. There are a source (S) and a receiver (R) at the surface. The wave travels vertically to a maximum altitude, where it is reflected, and it travels through a cross-wind field through all its trajectory. The top panel shows the original set-up with a native atmospheric grid. The middle panel shows the situation after reducing the problem to the along-track plane and discretising the (interpolated) cross-wind both in two directions. The bottom panel shows the problem after further simplifying by averaging the cross-winds on the along-track direction.
Many systems of interest have non-linear evolution and observation operators, i.e. m : R Nx → R Nx and h : R Nx → R Ny . This yields: The KF can still be implemented after linearising these operators. This is known as Extended Kalman Filter (EKF: see, e.g. Jazwinski [1970]) and its use is complicated since it involves large Jacobian matrices. An alternative is to use sample estimators for mean and covariance and to work with ensembles [Evensen 1994]. Hunt et al. [2007] nicely describe handling non-linear operators for this approach: Start with an ensemble of N e initial states, i.e. the matrix X 0 ∈ R Nx×Ne : The background ensemble at time t is found by applying the model to each member: We now drop the time index. The sample background meanx b ∈ R Nx is: A matrix of background perturbationsX b ∈ N x × N e is computed aŝ which relates to the (low-rank) sample covariance matrixP b as: The background ensemble in observation space Y b ∈ R Ny×Nx is obtained by applying the observation operator to each ensemble member of the background: The sample meanȳ b ∈ R Ny is computed as:ȳ This allows to compute the following expressions: and the computation of the ensemble-based gain is simplỹ