Convective gravity wave propagation and breaking in the stratosphere: comparison between WRF model simulations and lidar data

Abstract. In this work we perform numerical simulations of convective gravity waves (GWs), using the WRF (Weather Research and Forecasting) model. We first run an idealized, simplified and highly resolved simulation with model top at 80 km. Below 60 km of altitude, a vertical grid spacing smaller than 1 km is supposed to reliably resolve the effects of GW breaking. An eastward linear wind shear interacts with the GW field generated by a single convective thunderstorm. After 70 min of integration time, averaging within a radius of 300 km from the storm centre, results show that wave breaking in the upper stratosphere is largely dominated by saturation effects, driving an average drag force up to −41 m s−1 day−1. In the lower stratosphere, mean wave drag is positive and equal to 4.4 m s−1 day−1. In a second step, realistic WRF simulations are compared with lidar measurements from the NDACC network (Network for the Detection of Atmospheric Composition Changes) of gravity wave potential energy (Ep) over OHP (Haute-Provence Observatory, southern France). Using a vertical grid spacing smaller than 1 km below 50 km of altitude, WRF seems to reliably reproduce the effect of GW dynamics and capture qualitative aspects of wave momentum and energy propagation and transfer to background mean flow. Averaging within a radius of 120 km from the storm centre, the resulting drag force for the study case (2 h storm) is negative in the higher (−1 m s−1 day−1) and positive in the lower stratosphere (0.23 m s−1 day−1). Vertical structures of simulated potential energy profiles are found to be in good agreement with those measured by lidar. Ep is mostly conserved with altitude in August while, in October, Ep decreases in the upper stratosphere to grow again in the lower mesosphere. On the other hand, the magnitude of simulated wave energy is clearly underestimated with respect to lidar data by about 3–4 times.

L. Costantino et al.: Comparison between WRF model simulations and lidar data lation) events, which are not accurately modelled in the absence of gravity waves (Dunkerton, 1997;Hitchman and Leovy, 1998;Ray et al., 1998;Richter et al., 2010;Kawatani et al., 2010, Evan et al., 2012. Extratropics are also strongly affected by GW activity with an estimated monthly average negative forcing, increasing from a few up to about 100 m s −1 day −1 between 0.1 and 0.01 hPa, which drives the wind reversal around the mesopause (Lindzen, 1981;Holton and Alexander, 1999). This value has been confirmed, among others, by recent high-resolution simulations of Watanabe et al. (2008).
GWs may have a wide spectrum of horizontal scales that range from a few to hundreds of kilometres with periods from minutes to hours. General circulation models (GCMs), coupling troposphere and stratosphere for climate studies, have generally a coarse resolution in the stratosphere, between 2 and 5 • horizontally and 3 km vertically (Alexander et al., 2010). This resolution is fine enough to resolve Rossby waves but not small-scale GWs. Hence, the momentum forcing generated by unresolved waves is parametrized and constrained by large-scale observations of temperature and wind in the upper troposphere and middle atmosphere (Kim et al., 2003;Alexander et al., 2010). Traditionally, the GW drag (GWD) parametrizations used in climate and weather forecasting models aim to adjust the structure of winter jets and the horizontal temperature gradient. They were firstly based on the parametrization of orographic waves (Palmer et al., 1986), characterized by zero phase speed waves and generated by subgrid topography. In more recent times, as model tops increased up to the mesosphere, new GW schemes have become necessary to account for waves with non-zero phase speed, generated by other factors than topography. Our limited understanding of wave sources makes it difficult to validate the realism of such parametrizations. However, recent articles have tried to start constraining gravity-wave parametrizations with observations (e.g. Geller et al., 2013).

Purposes and strategy
To reduce uncertainties associated with GW parametrization in GCM, we need an improved knowledge of GW spectrum and its variability with altitude, wave sources, momentum and energy transfer, wave drag, wave breaking mechanisms and breaking heights. Regional mesoscale models, with horizontal resolutions that can reach a few hundreds of metres, are able to simulate small-scale GW activity. They may represent a valuable addition to direct ground-based (often limited in space and time) or space-based (often limited in resolution) observations, in order to analyse a large number of GW parameters. In this work, which is part of ARISE project (Atmospheric dynamics Research InfraStructure in Europe, http://arise-project.eu, an international collaborative infrastructure design study project funded by the FP7 European Commission), we make use of the mesoscale WRF (Weather Research and Forecasting) model (Skamarock et al., 2008; information online at www.wrf-model.org/index.php) to explicitly resolve wave motion (without any GW parametrization) and investigate GW propagation in the stratosphere and lower mesosphere.
Besides topography, it has been shown that deep convection is one of the most important GW sources in the stratosphere (i.e. Zhang et al., 2012). Here, we focus on GWs generated by deep convection, with the aim of quantifying horizontal momentum fluxes (HMFs) and wave drag forces above convective cells, as well as the amount of potential energy transported in the upper levels of the stratosphere. Alexander et al. (2010) showed that a minimum vertical resolution higher than 1 km is needed to reliably resolve GW activity throughout the middle atmosphere, together with a sufficiently high model top (near 1 Pa). According to these results, we run a bi-dimensional (2-D) and a highly resolved idealized simulation (idealized case), where a convective GW interacts with a stratospheric linear wind shear. Results are analysed and interpreted with respect to mediumand high-frequency GW linear theory. In a second step, we run more complex three-dimensional (3-D) simulations over southern France (real case). Realistic temperature and wind profiles, from ECMWF re-analysis data, are used as input values at model outer boundaries. Model results are then compared with co-located in situ observations of potential energy vertical profiles, measured by a lidar system at Haute-Provence Observatory (OHP) in southern France (43.93 • N, 5.71 • E).
The WRF model has been recently used to simulate real meteorological events and observe convective GW propagation in the stratosphere. These studies showed a good agreement between WRF simulations and observations. For instance, encouraging results come from Spiga et al. (2008), which used WRF to model inertial GWs (IGWs, with frequencies close to the inertial frequency) emitted above a convective cell in the lower stratosphere of Andes Cordillera region. With a resolution of 7 km and 500 m in horizontal and vertical grid spacing, respectively, their simulations allow the characterization of the sources of observed IGWs and the establishment of their link with vertical shears of horizontal wind. Comparing model results to ECMWF and NCEP-NCAR reanalysis, satellite and radio-soundings data radio, they clearly state good performance of the WRF model, which captured systematically the emitted IGW. With a coarser resolution (27 km of horizontal grid spacing), Kim and Chun (2010) simulated stratospheric gravity waves generated by a typhoon that moved in 2006 over the Korean Peninsula, showing a good agreement with both satellite observations and ECMWF analysis data. At larger scale (37 km of horizontal grid spacing), Evan et al. (2012) have been able to simulated convective GW in the ITCZ (Intertropical Convergence Zone), demonstrating and quantifying the role of GW forcing on QBO (Quasi-Biennial Oscillation).
In order to provide potential energy data to be compared with model simulations, Rayleigh lidar offers the unique pos-sibility to obtain temperature profiles in the upper stratosphere and lower mesosphere, with a high spatial and temporal resolution. Using lidar-based estimates of GW potential temperature, a number of field campaigns have confirmed the capability of lidar systems to capture main features of wave interaction with stratospheric mean flow. They showed a strong correlation between stratospheric wind intensity and E p . This correlation is found to be linear, with a correlation coefficient equal to 0.7 above southern France at 44 • N (Wilson et al., 1991) and Alaska at 65 • N (Thurairajah et al., 2010), and equal to 0.5 above Antarctica at 69 • S (Alexander et al., 2011). Alexander et al. (2011) analyse 839 h of temperature records during the autumns and winters of 2007 and 2008 and investigate the seasonal variability of GW (with vertical wavelength between 4 and 20 km and ground-based periods larger than 2 h). They find a peak in GW activity during winter up to an altitude of 40 km, above which the zonal wind starts to decrease and GWs propagate less efficiently. In autumn, GWs dissipate between 35 and 50 km but energy is conserved in the mesosphere.

Previous estimates of horizontal momentum fluxes
A number of previous works have already attempted to provide a quantification of GW momentum fluxes in the stratosphere, using high-resolution model simulations. However, theses studies are generally performed at global scale, with a horizontal resolution which is usually much coarser than that used in here. We provide hereafter some results, to give an order of magnitude of HMFs and drag forces already calculated in the stratosphere. However, the differences in resolution, domain size and simulation time window do not allow a direct comparison with HMF estimates provided in this study.
For instance, Sato et al. (2009) investigated the seasonal and inter-annual variations of GW in the stratosphere and mesosphere by a high-resolution global spectral model (the T213L256 middle atmosphere GCM developed for the KANTO project, see Watanabe et al., 2008), with a horizontal and vertical grid spacing corresponding to about 60 km (near the Equator) and 300 m (throughout the middle atmosphere), respectively. The authors state clearly that this horizontal resolution is insufficient to resolve very small gravity waves on the scale of 10 km. However, the vertical resolution is supposed to be sufficiently fine to resolve the majority of observed gravity waves with acceptable accuracy. HMFs exhibit an annual variation that is positive in summer and negative in winter (relative to each hemisphere), in both the lower stratosphere and mesosphere. In particular, in the winter Southern Hemisphere, they find large negative momentum fluxes near high mountains (Antarctic Peninsula, Andes, east coast of Australia) in good agreement with Plougonven et al. (2008) and Jiang et al. (2013), and distributed zonally along the eastward jet. In the summer Northern Hemisphere (NH), large positive momentum fluxes are observed only over Indian and African monsoon regions, probably generated by deep convection. Their results suggest different GW sources in winter and summer. Averaged zonally, the global monthly HMF varies between 2 (at the Equator, with peaks during the NH summer) and −4 × 10 −3 N m −2 (at tropics, with peaks during the winter of each hemisphere) at 100 hPa (lower stratosphere). At 0.03 hPa (upper mesosphere), HMF is close to zero near the Equator and oscillates between 0.12 and −0.2 × 10 −3 N m −2 (summer and winter of each hemisphere, respectively) at midlatitudes. These estimates are consistent with satellite data at 20-30 km, collected by EOS-Aura satellite and analysed by Alexander et al. (2008). During August 2006, they calculate absolute values of HMF decreasing from 1 (tropics) down to 0.2 × 10 −3 N m −2 (NH midlatitudes). Watanabe et al. (2008) used the KANTO model and found a positive GW forcing, in the extratropical NH during summertime, which increases with altitude from a few m s −1 day −1 (between approximately 3 and 0.01 hPa) up to 100 m s −1 day −1 (between 3 and 0.01 hPa), confirming the estimates of Lindzen (1981).
More recently, Geller et al. (2013) used CAM5 (horizontal resolution of 0.23 • latitude and 0.31 • longitude) and KANTO models (together with three other global models with much coarser resolution using GW parametrizations, not discussed in here) and compared them to satellite, balloon and radiosonde observations. Both models show a zonal mean of absolute momentum flux which is less than 2 × 10 −3 N m −2 at 45 • N, at 20 km of altitude, for July 2005, 2006 and 2007, in good agreement with results of Sato et al. (2009). However, in their conclusions, they stress how these two models under-resolved short-wavelength GW, associated with important momentum fluxes.
The importance of small-scale GW in the momentum budget of the stratosphere has been also stressed in a number of highly resolved and idealized experiments. For instance, Lane and Sharman (2006) used a three-dimensional model with horizontal and vertical grid spacing of 150 m. They showed that deep convective clouds generate GW of about 5-10 km of horizontal wavelength, followed by the occurrence of secondary smaller waves (2 km of horizontal wavelength) in the lower stratosphere (between 15 and 17 km of altitude) generated by the primary wave breakdown. Lane and Moncrieff (2008) extended the study up to 40 km. They observed that slower moving and short-scale tropospheric GWs (3-4 km of horizontal wavelength) make the strongest contribution to stratospheric vertical fluxes of horizontal momentum, even if the spectrum of GWs with the strongest power has horizontal wavelengths between 5 and 50 km.
All these results encourage the analysis of convective GW momentum and drag forces using highly resolved models, with horizontal grid spacing of at least 1 km, also when simulating real meteorological events.

The dispersion relation
In the following analysis, GWs will be treated as smallamplitude perturbation to some larger-scale horizontally uniform and steady background state. Considering only the xz plane, the dispersion relation for GW relates the frequency (ω) to the horizontal (k) and vertical (m) wavenumber. According to the derivation of Fritts and Alexander (2003), the simplified dispersion relation of the intrinsic frequency (i.e. the frequency in the reference frame of moving background atmospheric flow) for high-frequency (for which Coriolis force can be neglected) and hydrostatic (k 2 m 2 ) waves can be written as where U 0 is the zonal-mean wind speed (if we just consider the longitudinal plane), α is the angle between the vertical and the lines of constant phase, and N is the Brunt-Väisälä frequency. For a vertical propagating wave m and k are real and the intrinsic frequency is confined to the range 0 < ω < N.
Wave phase speed, c x , can be expressed as (Nappo, 2002) If we think at a bi-dimensional GW field propagating from its source, k is negative on the left side (with respect to GW source). At constant altitude, in the hypothesis that N, k and α do not vary significantly in space (along x) and time, and U 0 > 0, phase velocity is expected to be smaller in the western direction (k < 0) than eastern one (k > 0). Those altitudes where ω → 0, above which waves become unstable and break, are called critical levels. According to Eq. (1), if k and U are in the same direction, as U 0 increases isophase lines turn horizontally until α = π/2 (hence ω → 0 and m → ∞) and the vertical wavelength (λ m = 2π/m) goes to zero. Wave energy and momentum are transmitted to the mean flow and converted into small-scale turbulent motion. Another way for GWs to transfer energy and momentum to the mean flow is the so-called wave saturation mechanism, occurring when wave amplitude is too large with respect to vertical wavelength. Isophase lines become very steep and waves break. This process is a consequence of the tendency of wave disturbances to increase in amplitude with height (as density decreases), by a factor of e z/2H . According to the linear theory, there is no mechanism which could prevent GWs from growing indefinitely. However, this would lead to non-physical solutions (e.g. negative pressures). In nature the amplitude of disturbances is bounded and beyond a certain threshold waves break.

Drag force and energy
HMF is expressed as ρ 0 u w , where ρ 0 is the atmospheric mean density at a given altitude and u and w are the horizontal and vertical wind perturbation amplitudes. The quantity ρ 0 u w is conserved with height in the absence of wave dissipation. With increasing height, density changes or variations in wind shear or in static stability can cause GW to break so that momentum flux in no longer constant. Its vertical gradient is a measure of the force that dissipating waves exert on the mean flow. According to the zonal-average zonal wind momentum equation, the resulting drag force is directly related to the mean flow velocity as A change in the momentum flux with increasing altitude would result in a net acceleration (or deceleration) of the mean flow. The drag force of orographic waves is generally negative, slowing the wind speed. Convectively forced GWs can alternatively accelerate or decelerate the mean flow, dragging the wind in the direction of phase speed of breaking waves. The total energy per unit mass (energy density) E 0 (J kg −1 ) is a good proxy to measure the GW activity. It is defined as where T 0 is the background atmospheric temperature (average is spatial or temporal) and T the temperature perturbation. According to VanZandt (1985), energy repartition between E k and E p is supposed to be constant for medium-frequency waves (i.e. f 2 ω 2 N 2 ). In this case, GW total energy can be somewhat deduced by the potential energy alone, which can be derived from simple temperature measurements. On the other hand, E p /E k decreases as GW internal frequency decreases, with almost no temperature fluctuation in case of very long wave periods. In the absence of dissipation, wave energy (per unit mass) is supposed to increase proportionally to the square of the wave-induced disturbance amplitude, i.e. by a factor of e z/H , generally referred to as conservative growth rate.

Methods and experiments
We use the version 3.5.1 of WRF model. Governing equations and all parametrization schemes and numerical meth-ods of WRF are extensively explained in Skamarock et al. (2008), while those used in here are briefly described in Costantino and Heinrich (2014). Data from real case simulation are then compared with lidar data, from the ARISE campaign at OHP.

Idealized case
For this experiment we use a horizontal resolution of dx = 1 km for 2000 horizontal grid points and 451 vertical levels. Below 40 km of altitude, vertical spacing (dz) is fairly constant, varying between 90 and 130 m. Then, it increases almost linearly until model top (at approximately 80 km), where dz is equal to 7.8 km. In this way, vertical resolution is higher than 1 km below 60 km of altitude. At model lateral boundaries, we set open boundary conditions. The initial atmosphere is horizontally stratified and convection is triggered by a warm bubble (WB) with a maximum intensity of 3 K (switched off after the first temporal step) and horizontal and vertical radius of 4 and 1.5 km, respectively. It is placed in the middle of the domain (x = 1000 km) at an altitude of z = 1.5 km. To reduce wave reflection at model top, in the last 10 km we put a Rayleigh absorbing layer with damping coefficient of 0.02 s −1 . For time integration, we use the third-order Runge-Kutta scheme with a time step of 2 s. Vertical and horizontal eddy diffusion coefficients are set to zero (i.e. no subgrid turbulence). Cloud dynamics is supposed to be fully resolved by motion equations, and no cumulus parametrization scheme is used. Initial meteorological parameters (humidity and potential temperature vertical profiles) are only a function of altitude and are derived from real case output (see next paragraph) over the Mediterranean Sea off the coast of southern France, very close to a precipitation field (11:00 on 21 October 2012). In that way, we try to be as close as possible to realistic extratropical NH conditions, during autumnal rain events. The initial background wind is a simplified profile of real case stratospheric wind shear. It is set to zero below 16 km, increasing linearly up to about 75 m s −1 at 65 km of altitude.

Real case
We perform two real case simulations for 21-25 August and 19-24 October 2012. We use a two-way nesting for three concentric domains over Europe and northern Africa (mother domain or grid 1), France and western Mediterranean Sea (grid 2), and southern France (grid 3), with horizontal resolutions of 27, 9 and 3 km, respectively. We use 131 vertical levels for October and only 101 for August, as meteorological conditions (strong winds over high mountains) were unfavourable for a finer resolution. In the 131-level run, the vertical resolution increases linearly from 50 to 300 m in the troposphere and remains almost constant up to 45 km of altitude. Above, dz increases linearly from 300 m to 4 km until model top (at 68 km and 7.3 Pa). In this way, vertical resolution is higher than 1 km below 50 km of altitude. The reference temperature in the stratosphere is 220 K. The 101-level simulation has very similar characteristics. The main difference with respect to the 131-level case is the lower vertical resolution of the last few levels, with dz increasing from 1 km (at 50 km) up to 5.7 km (at 66 km).
Kain-Fritsch convective parametrization scheme is applied only to grid 1 and 2. The relatively small resolution of grid 3 is supposed to resolve cloud dynamics and cumulus parametrization is not used, as suggested by Skamarock et al. (2008). For cloud microphysics, we use the Ferrier scheme. The boundary layer scheme is that of the Yonsei University (YSU).
The principal (meteorological) time step of the third-order Runge-Kutta integration scheme is equal to 30, 10 and 2.5 s (for the three different grids respectively), while the secondary time step (resolving acoustic waves) is 256 times smaller. The model is fully non-hydrostatic. Coriolis force acts only on wind perturbations. Turbulent eddy coefficients are calculated using the horizontal Smagorinsky first-order closure.
Note that the three domains are not nudged. ECMWF reanalysis data are only used to provide realistic meteorological conditions (horizontal wind components, temperature and specific humidity) at mother grid boundaries every 3 h.
At model top, the w-Rayleigh layer depth is 10 km, with a relatively high damping coefficient equal to 0.2 s −1 . Topography has been smoothed in both WRF and WPS (WRF Preprocessing System).

Instrument description
The lidar instruments are powerful tools for the study of atmospheric perturbations. They produce accurate observations with high temporal and spatial resolution, well adapted for studying atmospheric gravity waves (Chanin and Hauchecorne, 1981). Gravity wave activity has been extensively analysed using lidar throughout the middle atmosphere in several studies (Fritts and Alexander, 2003, and references therein).
Rayleigh lidar provides vertical profiles of molecular density and temperature when the atmosphere is free of aerosols (Rayleigh scattering above 30 km) from about 30 to 90 km depending on the signal-to-noise ratio (Hauchecorne and Chanin, 1980). The OHP lidar is composed of a frequencydoubled Nd:YAG laser emitting at 532 nm with a repetition rate of 50 Hz and a collector surface area composed by a mosaic of four mirrors with a diameter of 50 cm corresponding to a surface of 0.8 m 2 . Lidar measurements have been performed continuously at OHP since late 1978. In early years, the vertical resolution was 0.3 km, and it has been improved to 0.075 km since the mid-1990s. The temporal resolution is about 2 min 40 s. We only used night-time profiles during clear-sky conditions above OHP. During night-time the background noise decreases considerably (for more details about the lidar instrument and technique, see Keckhut et al., 1993). The number of observations used in this work is 14 profiles for August and 13 profiles for October 2012.

Variance method: a brief description
In order to have access to perturbations with short time and vertical scales, at least in a statistical sense, we analyse raw lidar signals with a variance method (Hauchecorne et al., 1994;Mzé et al., 2014). This method is based on the computation of the signal perturbations over short time and vertical intervals ( t, z) and on the summation of the square of these perturbations over a large number of elementary intervals ( T = tN t , Z = zN z ), which give an estimation of their variance. It allows extracting root-mean-square mean amplitude of small-scale perturbations that are not detectable on single profiles. The observed variance of the signal is defined as the sum of instrumental and atmospheric variances: The instrumental variance is estimated assuming a Poisson noise distribution for the photon counting signal. Then, the atmospheric variance will provide an estimation of the GW activity in the middle atmosphere.
The estimation of the variance with a given thickness Z = zN z of the layers is equivalent to estimating the power spectral density of atmospheric fluctuations in bandpass filter with characteristics related to Z. This method is equivalent to an estimation of the variance using a broad band-pass filter centred at λ v ∼ 2.4 Z. Between 30 and 50 km Z =1.5 km and N z = 20, while Z = 3 km and N z = 40 above. Lidar vertical resolution ( z) is equal to 0.075 km. T = tN t ∼ 26.7 min, with t = 160 s and N t = 10.
The variance method is computed for each night between 30 and 85 km of altitude and expressed in potential energy per unit mass (J kg −1 ) using the Brunt-Väisälä frequency from the mean lidar temperature profile, in order to characterize gravity wave activity.
The variance method is relatively simple but has the advantage of being robust, fast and using raw data. It is independent of data processing errors. More details about the variance method are presented in a recent study by Mzé et al. (2014).

GW propagation
A number of sensitivity studies have shown that the sponge layer has an evident strong damping effect on the zonal flow forcing. For that reason, in the following analysis we con-sider only the region below 68 km of altitude (indicated by a dashed line in figures).
Under the effect of the initial instability generated by the WB at t = 0, the first cloud forms after about 4 min integration time, and precipitation after 8 min. In good agreement with the results of Costantino and Heinrich (2014), the cloud grows up quickly and reaches the tropopause (12.5 km of altitude) after about 20 min. At this point, a large spectrum of gravity waves is generated close to the tropopause and propagates downward (in the troposphere) and upward (in the stratosphere) direction. The initial wind is a simplified profile close to that observed near the OHP in the real simulation, on 21 October 2012, at 00:12 UTC. It is zero up to 17.5 km (in the real simulation it is slightly negative) and then it increases linearly up to 80 m s −1 at 65 km (in the real simulation the wind shear is not linear, but the intensity is similar). Figure 1 shows the horizontal cross section of wind vertical velocity, w, which is a good proxy to observe gravity waves. Green contour lines represent absolute values of the ratio between potential temperature anomaly and its horizontal mean value, by steps of 0.02. With the term anomaly we define the difference between the local value of a physical quantity and the horizontal average. This parameter is an indicator of atmospheric buoyancy, which is a source of instability and turbulence. Figure 1 shows a continuous generation and propagation of gravity wave packets from 40 to 90 min of integration time, by steps of 10 min. In the stratosphere, energy is transported upward and eastward on the right side of the domain, westward on the left side. This can be seen following the temporal displacement of vertical velocity intensity peaks, which is reliably representative of GW group velocity.
There is a strong asymmetry in wave propagation between eastern and western part, which is driven by the eastward wind shear. According to Eq. (1), a change in the horizontal wind speed has a direct effect on α (the angle between the vertical and the isophases), which decreases (increases) with altitude on the left (right) side of the storm.
On the left side (upwind) α remains small but positive ,and waves can propagate vertically in the stratosphere (until α is zero and waves become evanescent). Their amplitude, increasing with altitude, may reach the threshold value beyond which wave saturation can occur. On the right side (downwind), α increases with altitude (mostly visible in the lower stratosphere) and wave intrinsic frequency decreases. According to Eq. (1), when U 0 is strong enough so that ω → 0, waves break and can no longer penetrate into the upper stratosphere. This process is generally referred to as wind filtering.
Horizontal phase velocity of wave packets at a given altitude can be analysed in more detail looking at Fig. 2, where the vertical wind speed is shown as a function of time (y axis) and space (x axis). Slopes of isophases represent the horizontal velocity of wave packet phases. As a reference, dark green Figure 1. Vertical cross section of wind vertical velocity, w (m s −1 ), cloud water mixing ratio (g kg −1 ) and rain water mixing ratio (g kg −1 ) at t = 40, 50, 60, 70, 80, and 90 min (from top to bottom). Arrow's length is proportional to the horizontal component of background wind velocity (1 km on the x axis is equal to 4 m s −1 ). Green contour lines show absolute values of the ratio between potential temperature anomaly and its horizontally averaged value. lines overplotted on figures show horizontal velocity slopes relative to 16.7 (dotted), 33.3 (solid) and 66.7 (dashed) m s −1 (i.e. 60, 120, 240 km h −1 , respectively). On the left side of the domain, horizontal phase speed of GW generated between 40 and 60 min of integration time, and within a radius of 100 km from storm centre, is close to 16.7 m s −1 . Horizontal wavelengths are equal to about 2-3 km, at 40 km of altitude, and 4 km at 65 km. This may give us a rough estimate of wave periods of a few minutes. Oldest waves, generated in the first 20 min of GW activity (20-40 min of integration time), have much stronger phase velocities (up to approximately 66.7 m s −1 ) and longer horizontal wavelengths, up to 10 km. This is in good agreement with the results of Lane and Sharman (2006) and the order of magnitude of convective GW wavelengths they found, between 2 and 10 km.
Top image of Fig. 2 (40 km) shows how phase velocity of wave packets may vary with time. For instance, at t = 60 (or t = 80) min and 900 < x < 950 km, isophase slopes rapidly increase, turning almost vertically (phase speed deceleration), and decrease again little afterwards (acceleration). As expected, phase velocity appears higher downwind than upwind, and also more constant with time. Comparing top and bottom images of Fig. 2, its seems that phase velocity is likely to be higher, on average, in the upper levels of the stratosphere. Figure 3 shows the horizontal momentum flux, calculated averaging the u w product over an area of 300 km on the left (blue line) and right side (red line) of the storm.

Momentum flux
The averaged value over the whole 600 km region (left and right side together) is reported in green. Before the generation of convective GW, HMF is almost zero in the stratosphere, while after t = 20 min, GW dynamics induces a vertical transport of horizontal momentum. At t = 30 min, HMF is negative upwind, where u w is generally negative (air parcels oscillate along a slant path in a upward-westward direction), and positive downwind. Stratospheric peak values of HMF are equal to −73 × 10 −3 and 76 × 10 −3 N m −2 at 14.6 km of altitude. With increasing altitude HMF tends rapidly to zero, with absolute values smaller than 2 × 10 −3 N m −2 beyond 40 km on the left side and 30 km on the right side. Upwind gravity waves, no subjected to wind filtering, penetrate deeper in the stratosphere. Considering the whole region, the balance between negative and positive values is slightly negative, with a peak value in the stratosphere of −6 × 10 −3 N m −2 (at 23.5 km).
With increasing time, the maximum altitude of absolute HMF values larger than 2 × 10 −3 N m −2 increases up to 49 (left side) and 38 km (right side). At t =50 min, it reaches 60 (left side) and 37 km (right side). Considering the whole region, the overall HMF is positive between 13.5 and 22 km of altitude (with a maximum of 12 × 10 −3 N m −2 , at z = 14 km) and negative above (with a peak of −12 × 10 −3 N m −2 , at z = 29 km). At t = 70 min, HMF is positive between 17.5 and 28.5 km of altitude (with a maximum of 17 × 10 −3 N m −2 , at z = 23 km) and negative above (with a peak of −8 × 10 −3 N m −2 , at z = 43 km).
In conclusion, the upward transport of HMF is very efficient, with an average vertical speed of about 30 km h −1 . In 40 min (from t = 30 to 70 min) the peak of negative HMF moves by 20 km (from z = 23 to 43 km), as well as the highest altitude with HMF smaller than −2 × 10 −3 N m −2 (from z = 40 to 60 km).
Note that we are using here the hypothesis that storm anvil, and hence GW source, is punctual and centred in the middle of the domain. This is not completely exact, as we can see in Fig. 1, showing that storm edges move outward up to 30 km away from domain centre (t = 70 min). Cloud edges are the most convective parts of a storm and GW are mostly generated at storm borders. Then, in a bi-dimensional simulation, we have actually two GW sources, very close and symmetric. In the region located just above cloud top (e.g. at t = 70 min, Fig. 1) a complex interaction between the downwind wave (propagating rightward from the left edge of the storm) and the upwind wave (propagating leftward from the right edge of the storm) occurs. Looking closely at the time development of this interaction field, it really seems that it does not propagate beyond the outer edges of the cloud. In particular the downwind wave is not advected eastward beyond the right edge, remaining bounded in the region within ±30 km, with a zero phase speed (this is maybe due to the interference of the two waves that are symmetric, with very similar characteristics and intensities). To test the consistency of our approximation (i.e. to consider the storm as punctual), we calculated HMF eliminating the area within ±30 km. Above 20 km of altitude, results are very close to those shown in Fig. 3, where the area within ±30 km is considered. We believe that our hypothesis, always remaining an approximation, does not alter considerably the results of our analysis. For real case simulations, this approximation is supposed to work even better as convective cloud anvils are smaller.
Note also that another factor (other than wave breaking) that may contribute to the decrease of eastward momentum in Fig. 3 is the fact that eastward-propagating waves leave the ±300 km subdomain (where momentum fluxes are diagnosed) by its lateral side. To verify this point, we have calculated the HMF in almost the entire domain (±900 km). Above 20 km of altitude, these results are very close to those obtained averaging over a subdomain of 300 km. In particular, qualitative features of HMF vertical profiles (e.g. such altitudes where HMF decreases or increases) are very similar for the two cases, indicating that a ±300 km averaging window is large enough not to lose essential information. Mean HMFs, averaged over the whole 600 km subregion, have the same order of magnitude of those calculated by a highly resolved global model in the lower stratosphere of midlatitude NH, of a few 10 −3 N m −2 (see Sect. 1.2). In our case, values are stronger, probably because measurements are only performed during storm occurrence. In addition, the smaller grid spacing of our simulation may resolve small GW scales associated with important HMFs that are neglected in global models. Figure 4 shows the horizontal drag force calculated at t = 70 min, on the left (blue line, left image) and right side (red line, middle image) of the domain, and the overall average value of left and right side together (green). As expected, downwind wave breaking leads to a strong negative peak of instantaneous drag force, down to −140 m s −1 day −1 at 52.5 km of altitude. The thinner black line shows the mean horizontal wind anomaly with respect to the left half of the domain. U-wind anomaly has a peak of −0.75 m s −1 at 52 km, coincident with that of drag force, and shows a similar vertical profile to GW drag distribution. Upwind wave breaking is significantly weaker and the resulting drag force (red line) is positive but very small, with two peaks in the lower and two in the upper stratosphere, at 18.5 km (13.6 m s −1 day −1 ), 25 km (11.9 m s −1 day −1 ), 54.3 km (6.8 m s −1 day −1 ) and 61.4 km (8.5 m s −1 day −1 ). The resulting acceleration of the mean wind is almost irrelevant, oscillating between positive and negative values. Near the tropopause, approximately below 15 km, the large drag force values and wind anomalies on both left and right sides are probably due to outward air flux above cloud top, driven by storm convective dynamics.

Drag force
Considering the whole area, the mean stratospheric drag force is slightly positive between 20 and 40 km (vertical average of 4.37 m s −1 day −1 ) and strongly negative above (ver-

Energy
Potential energy is calculated according to Eq. (6), averaging spatially over the whole domain (left and right side together). Figure 5 shows vertical profiles of potential (blue), kinetic (red) and total energy (black) at t = 30, 50 and 70 min of integration time. Green dashed line represents conservative growth rate (Sect. 2.3). As expected, the development of the storm and the consequent excitation and propagation of GWs increase enormously the wave energy in the stratosphere. At t = 30 min, mean total energy between 20 and 60 km is very small and equal on average to 1 J kg −1 . Above the tropopause, its vertical profile is almost constant with height, a sign of a small GW penetration in the upper levels. At t = 50 min, stratospheric total energy is 10 times larger, on average, than 20 min before (11 J kg −1 ) and attains a maximum value of 20 J kg −1 at 55 km. At t = 70 min, the growth rate is almost exponential between 30 and 50 km. Mean total energy of the stratosphere reaches 30 J kg −1 , with a strong peak of 68 J kg −1 at z = 50 km. Above this altitude, the sudden decrease in E p suggests the presence of wave breaking by saturation. This mechanism is supposed to increase atmospheric instability and may explain the further increase in kinetic energy for z > 50 km. The close profile and mean values of potential and kinetic energy confirm, to a certain extent, the hypothesis of constant energy repartition between E k and E p (Sect. 2.3) in particular below 40 km.

Meteorology and GW dynamics
Similarly to the idealized case, in the following paragraphs we focus on the study of stratospheric dynamics up to about 58 km of altitude, where the damping layer starts. This altitude is marked in figures with a dashed line.
The case study is the 19-24 October experiment. In particular, on 21 and 22 October, a very strong rain thunderstorm occurs in the western Mediterranean Basin. From a depression over eastern Spain, cold air converges with warm and humid air that flows from northern Africa toward southern and eastern France and Germany. Over France, a cold front forms and moves eastward from Pyrenees (Spain-France border) to south-east. Figure 6 shows infrared data acquired by the IR channel (at 10.8 µm) of the geostationary meteorological satellite METEOSAT9 and provided by EUMET-SAT (European Organization for the Exploitation of Meteorological Satellites). From left to right, images are relative to 21 October, at 14:00 and 21:00 UTC, and 22 October at 04:00 UTC. The IR 10.8 µm channel provides information on cloud top temperature. Colder cloud tops, generally associated with more vertically developed clouds, are whiter in the colour scale. Figure 6 shows a large mesoscale cloud system which is advected eastward over the western Mediterranean Basin. The bright and white spot in southern France at 21:00 and 04:00 UTC suggests the presence of strongly convective clouds probably associated with severe weather conditions over the OHP region (indicated in figure by an orange square).
Temporal and spatial coincidences between real and modelled meteorology are fairly good. According to WRF simulation, in the morning of 21 October (from 09:00 UTC) a heavy rainstorm appears in grid 3, from the south-east boundary. It passes twice over the OHP station from 16:00 to 18:00 UTC (21 October) and from 20:00 (21 October) to 04:00 UTC (22 October), with the strongest rain event be-tween 22:00 and 00:00. Figure 7 shows the horizontal cross section of vertical velocity (grey colour scale) at 10 and 40 km of altitude, for 21 October, at 23:00 UTC. Contour lines represent the rain water mixing ratio (qrain), by steps of 1 g kg −1 . The maximum of the total column (from the ground to top of the atmosphere) qrain value (at 3 km resolution) is reported under the image, and it is equal to 44.4 g kg −1 . The position of OHP station is shown in green.
In good agreement with satellite observations, Fig. 7 shows a strong rain storm occurring over the OHP region. The severe convective dynamics (clearly visible in Fig. 8) perturbs significantly the vertical velocity field. While at low altitude levels (left image), the strongest w-speed perturbations are mostly coincident with the highest peak of rain (central part of the domain) and the highest mountain (as the Alps, in the north-east of the domain), at higher altitudes w field is more coherent and shows a large spectrum of gravity waves that propagate in the stratosphere to very long distances. The yellow square represents the grid box (referred to as virtual station, VS) where the product of vertical velocity and qrain is maximum, i.e. the place where storm convective dynamics is supposed to have the largest impact on the atmosphere. At 40 km of altitude, the concentric rings of w crests seem to indicate that the most active GW field originates from the virtual station and propagate upwind, in the westward direction. Downwind (eastward), GW propagation appears much less efficient, probably because of the wind filtering by wave breaking.
GW dynamics is shown in deeper detail in Fig. 8, representing a vertical cross section of zonal (top image) and meridional (bottom image) vertical velocity. The blue-red colour scale is for the rain water mixing ratio, while the grey scale is for the cloud mixing ratio. The vertical black line indicates the location of VS. Wind vectors are overplotted and arrow's length is proportional to the horizontal wind intensity (0.1 • of latitude or longitude is equal to 10 m s −1 ). Deep convection occurs when cloud top reaches the tropopause. This is almost the case above the VS, where the tropopause is at 12.5 km and cloud top attains 10 km of altitude. Above cloud top, the ascending flow has a strong positive vertical speed and interacts with the tropopause, from where a large GW field originates. This mechanism of wave generation is very close to that observed in the idealized case. In the top image (zonal cross section), GWs propagate upward and outward, with respect to the storm. According to linear theory, wave propagation is much more efficient upwind, where the strong zonal wind (up to 80 m s −1 at 70 km) turns the isophases. The wave packets reach the left border of the domain, 300 km westward (the distance in kilometres from the left boundary is indicated on the upper x axis). Downwind, isophases turn horizontal with increasing altitude up to 45 km (critical level), where vertical wavelengths go to zero and GWs disappear. Similarly, meridional wave propagation (bottom image) is much more efficient upwind (northward) than downwind. On the other hand, meridional wave activity is weaker Figure 5. Spatial average of potential (blue), kinetic (red) and total (black) energy (J kg −1 ), calculated using Eq. (6), within a radius of 300 km from storm centre, at t = 30, 50 and 70 min of integration time. Vertically averaged values of total, potential and kinetic mean energy are reported in figure, in J kg −1 . Green dashed line represents conservative growth rate. The colour scale, from black to white, is proportional to cloud top temperature. Colder cloud tops (generally associated with more vertically developed clouds) are brighter. The OHP station (5.7, 43.9 • ) in southern France is indicated by a red-yellow square. than a zonal one as wind intensity is lower and close to zero, with punctual reversals at 30, 40 and 55 km. On the right side of the image, above 44 • N (northern part of the domain), the presence of high mountains (Alps) together with a strong northward low-tropospheric wind generates intense orographic waves. The wind reversal above the tropopause, however, prevents them from propagating further southward and upward into the lower stratosphere.

Momentum flux, drag force and energy
We quantify the drag force exerted on the mean flow by the thunderstorm that occurred during the night of 21 October, from 22:00 to 00:00 UTC. During this time period, the storm attains its maximum intensity and it is approximately located above the OHP. For each WRF 5 min output, HMF is calculated averaging horizontally the u and w product over a square of ±40 grid points (i.e. ±120 km in the S-N and E-W directions) from OHP position. From the vertical derivative of HMF we obtain the instantaneous drag force, which is then averaged temporally over the 2 h time period. Error bars represent the confidence level of mean values. They are calculated as σ/(n − 2) 1/2 , where n is the number of instantaneous measurements within each altitude bin and σ is their standard deviation. Spatially and temporally averaged drag force vertical profiles for the left side (blue), right side (red) and the whole square region (green) are shown in Fig. 9.
Above 16 km, stratospheric zonal wind (solid black line) is directed eastward and increases with altitude up to 65 m s −1 at 60 km. Meridional wind (dashed black line) is particularly weak and changes direction several times with increasing altitude. This is a typical meteorological condition in the Northern Hemisphere during autumn and winter. In good agreement with the idealized case, upwind drag force is small in the lower stratosphere, where no wave saturation effect is expected. A double negative peak equal to −2.31 and −3.60 m s −1 day −1 occurs at 41.2 and 45.1 km of altitude. Then the deceleration approaches to zero at 50 km of altitude, increases beyond this altitude and attains −6.5 m s −1 day −1 at 56.1 km, which is the highest altitude level not directly affected by the damping. Sensitivity stud- Arrow's length is proportional to the horizontal wind intensity (0.1 • of latitude or longitude is equal to 10 m s −1 ). Contour lines represent the total column rain water mixing ratio, qrain (g kg −1 ), by steps of 1 g kg −1 . Figure 8. Zonal (top image) and meridional (bottom image) vertical cross section of wind vertical velocity, w (m s −1 ), cloud water mixing ratio (g kg −1 ) and rain water mixing ratio (g kg −1 ). Arrow's length is proportional to the horizontal wind intensity (0.1 • of latitude or longitude is equal to 10 m s −1 ).
ies with different sponge layer depths show that the strongly negative drag force values in the last 10 km layer (down to −51.7 m s −1 day −1 at 65.3 km) should be considered as an artefact due to the wind deceleration imposed by the sponge layer.
Downwind forcing shows two peaks in the lower stratosphere, at 25.5 and 32.3 km, of 1.84 and 1.74 m s −1 day −1 , and a stronger peak in the stratosphere of 4.3 m s −1 day −1 , at 42.0 km. Also in this case, there is a good qualitative agreement with the idealized case, with positive peaks both in the lower part and the higher part of the stratosphere.
Averaged over the whole region (green line), temporal mean of drag force vertical profile has two positive peaks of 1.38 and 1.39 m s −1 day −1 (at 31.7 and 43 km) and two negative peaks of −1.87 and −4.53 m s −1 day −1 (at 45.6 and 56 km). Averaged vertically, drag force is positive in the lower stratosphere between 20 and 40 km (0.23 m s −1 day −1 ) and negative in the upper layers between 40 and 58 km (−1.00 m s −1 day −1 ).
These values are more than 10 times smaller in magnitude than those observed during the idealized storm up to 17 × 10 −3 N m −2 (23 km) and −8 × 10 −3 N m −2 (43 km), while vertical distribution of positive and negative maxima is in very good agreement. On the other hand, they are consistent with previous model results presented in Sect. 1.2.
If we look at the instantaneous energy (Fig. 10, right image), we arrive to similar conclusion. The intensity of instantaneous real case wave energy is, however, less underestimated than HMF. With respect to the idealized case at t = 70 min, real case energies (E p and E k ) on 21 October at 23:00 UTC are between 3 and 4 times smaller. Total wave energy is increased by about 58 % with respect to no rain conditions (21 October at 08:00 UTC), from 6.67 (E p = 3.55, E k = 3.12) to 10.56 (E p = 4.83, E k = 5.73) J kg −1 .

Comparison with lidar data
Here we compare the potential energy calculated by WRF for the two study cases of 21-25 August and 19-24 October 2012, with lidar measurements collected during the months of August and October 2012. For model results, we average over the whole time period the instantaneous potential energy relative to the single OHP grid box. Lidar data are collected only in case of clear sky in a narrow time window of about 3 h.
The variance method, computed for each night, provides one average profile per day from data collected approximatively in the 19:00-22:00 UTC time window. The variance is computed with N z = 20, between 30 and 50 km, and N z = 40 between 50 and 85 km. N z defines the average wavelength selected by the band-pass vertical filter used. N z = 20 corresponds to a band-pass filter centred at about 3.6 km (with a spectral interval between 2.4 and 5.8 km at half maximum), while N z = 40 corresponds to a band-pass filter centred at about 7.1 km (with a spectral interval between 5.1 and 11.3 km at half maximum). Single lidar profiles (provided every 160 s) have been integrated over ∼ 26.7 min (N t = 10) and then averaged over the 3 h time window. Stratospheric convective waves are deep (wavelengths of ∼ 10km). This corresponds to the part of the spectrum of gravity waves we intend to capture with lidar, using a spectral window centred at 7.1 km in the upper stratosphere, consistent with the order of magnitude of GW wavelength reproduced in our simulations. The integration time used limits the shortest period which can be measured (∼ 52 min).
To reduce statistical uncertainties due to the small number of available lidar data for the 21-25 August and 19-24 October time period (just 3 and 2 daily profiles, respectively), we consider all measurements collected during August and October for a total of 14 and 13 daily profiles, respectively. Lidar data are then supposed to describe, to a certain extent, the average effect of monthly atmospheric dynamics on GW energy. Error bars indicate the confidence level, with respect to the temporal variability, in both WRF and lidar data. Figure 11 shows in red WRF mean E p profiles together with lidar monthly averages (blue), for August (left) and October (right). The green lines represent the conservative growth rate, while the black horizontal line at 58 km of altitude indicates the beginning of the sponge layer. During August, lidar data show that E p mostly follows a conservative increases from 34 to 62 km of altitude, with a punctual but sensible energy drop between 52 and 54 km. In October, lidar measurements show that vertical energy transport is very efficient and conservative only between 36 and 44 km. Above this altitude there is a strong departure of E p from the conservative growth rate up to 52 km, where E p starts again to increase with increasing altitude but much less than in August. In the last four altitude levels, above 57 km, the average E p attains 15.4 J kg −1 in August and 9.7 J kg −1 in October.
WRF seems to capture the main E p feature revealed by the lidar system. In August, model results show an exponential increase of energy with altitude above 30 km in August, while in October there is a clear energy loss between 44 and 50 km, somewhat coincident with that seen in the lidar profile.
In Fig. 12, we present the scatterplot of WRF and lidar data for each month. In August (left image), the reduced chi square (i.e. divided by the degrees of freedom, df) is very close to 1. This value indicates a good agreement between the two data sets, with respect to the error variance. In October, where the energy vertical profile is much more variable with altitude, the reduced chi square is slightly higher but still close to the unity and equal to 1.67. However, despite a relatively good linear correlation between the two data set, WRF seems to underestimate systematically lidar values by a factor of 3 (October) and 4 (August).  (blue), right side (red) and whole domain (green), within a square of ±40 grid boxes (i.e. ±120 km on both x and y direction) from OHP. Right panel: spatial average of potential (blue), kinetic (red) and total (black) energy (J kg −1 ); 21 October 2012, at 23:00 UTC of (coincident with the peak of thunderstorm intensity).
Vertically averaged values of total, potential and kinetic mean energy are reported in J kg −1 . Green dashed line represents conservative growth rate.
OHP is located very close to Alps and mountain waves can have a strong impact on the atmospheric energy budget measured by lidar, which works in clear-sky conditions. At the same time, WRF is capable of resolving mountain waves and their energy is supposed to be fully captured, at least in the inner domain. However, U and V wind profiles shown in Fig. 9 indicate that over the study region wind inversion at the tropopause prevents the largest part of mountain waves from propagating upward. Hence, orographic GWs are not supposed to contribute consistently to the stratospheric energy budget near OHP. Outer domains, however, have a coarser resolution. Mountain and convective GWs occurring outside the inner borders are only partially resolved and their energy can not totally propagate into the inner domain. During the study period, large weather perturbations occur all around Europe (Fig. 6). Over the OHP, the lack of energy contribution from convective (and eventually orographic) GWs generated outside inner boundaries can be a leading factor of WRF energy underestimation with respect to lidar data.

Summary and conclusions
In this work we perform idealized and real case simulations of gravity waves generated by thunderstorms during the summer and autumn of midlatitude Northern Hemisphere. With respect to real case, the idealized model uses a simplified framework (e.g. flat orography, constant wind shear) that allows higher model top and higher horizontal (1 km) and vertical grid spacing (less than 1 km below 60 km of altitude). In the inner grid, the real case uses a horizontal resolution of 3 km and a vertical one of less than 1 km below 50 km. Previous studies (e.g. Alexander et al., 2010) suggest that a vertical grid spacing smaller than 1 km is fine enough to resolve the effect of GW breaking. Here, we test the capability of WRF to reproduce atmospheric dynamics up to 60-70 km of altitude, during deep convective rain events. We first analyse the idealized experiment and then we compare real case simulations with real energy data from lidar measurements collected during the ARISE campaign over the OHP station in southern France.
In the idealized case, convection is triggered by a warm bubble of 3 K within a very stable environment. In the real case simulation, we study two storm events during the periods of 21-25 August and 19-24 October 2012.
Even if background thermodynamic conditions are similar, idealized and real case experiments are intrinsically different and not directly comparable. First of all, that is because convective GW sources (and their magnitude) are not the same. In addition, 2-D numerical experiments are supposed to overestimate the energy of the whole meteorological system, because of lack of energy loss in the third dimension. Moreover, it has been shown that wave breaking itself is expected to be a highly three-dimensional process (Andreassen et al., 1994), even if qualitative aspects (momentum flux distribution and drag force) are fully captured in 2-D simulations. Finally real case has a much coarser resolution that can be a serious limit to reproduce the same amount of GW momentum and energy in the stratosphere.
We observe that in both idealized and real experiments deep convection is a very efficient source of small-scale GWs that propagate from the tropopause up to 60 km. Energy and horizontal momentum fluxes are transported from below to the high stratosphere and lower mesosphere.
For what concerns HMF, in the idealized simulation the maximum attains −8 × 10 −3 N m −2 at 52 km of altitude, which is 10 times larger than that observed in the high stratosphere on 21 October 2012, equal to −0.8 × 10 −3 N m −2 .
In the idealized case, drag force vertical profile is negative in the upper stratosphere (with an average value of −40.9 m s −1 day −1 ) and slightly positive in the lower  Scatterplot between real case WRF simulation (x axis) and lidar data (y axis) of potential energy per unit mass (J kg −1 ) for August (left) and October (right) 2012. Dotted line represents the best line considering errors on both x and y. The equation of best-line fit is reported, together with the reduced chi-square value (i.e. divided by the degrees of freedom, df). Error bars as in Fig. 12. stratosphere (4.4 m s −1 day −1 ). The peak value attains −68 m s −1 day −1 at 53 km of altitude consistent with expectations (Lindzen, 1981;Holton and Alexander, 1999). On the other hand, if we consider the real meteorological case (where a long transversal squall line with multiple centres passes above OHP, moving eastward with time) values are strongly different in magnitude. Note that because of storm motion, upwind and downwind momentum fluxes (calculated with respect to OHP) may result somewhat mixed so that spatial average of GW forcing can be underestimated. During the strong rain events of 21 October (22:00-00:00 UTC), the real case mean drag force (within a radius of 140 km from OHP) has two positive peaks of 1.38 and 1.39 m s −1 day −1 (at altitudes of 31.7 and 43 km) and negative values of −1.87 and −4.53 m s −1 day −1 (at altitudes of 45.6 and 56 km). On average, the forcing is positive in the lower stratosphere (0.23 m s −1 day −1 ) and negative in the upper layers (−1.00 m s −1 day −1 ). The presence of a strong damping layer above 58 km in real case experiments does not allow reliable drag force estimates at mesospheric altitudes.
In conclusion, the magnitude of real case HMF and wave drag force over OHP are consistent with previous NH midlatitude monthly average from Alexander et al. (2008), Sato et al. (2009), Watanabe et al. (2008 and Geller et al. (2013). Our results show slightly higher HMF. This can be due to the fact that we study punctual storm events in a narrow time window, but also due to higher horizontal resolution of our simulations that can resolve smaller-scale GW associated with stronger fluxes.
For what concerns energy, stratospheric E p profile of idealized case indicates a conservative growth rate up to 50-55 km, where wave saturation seems to be the strongest mechanism in decreasing wave energy. On the other hand, real case experiments show a strong dependence of E p on the time period, very consistent with in situ observations. Lidar vertical profiles of GW potential energy (per unit mass) show that in August the monthly averaged E p is almost completely conserved in the stratosphere, with an exponential increase with altitude very close to the conservative growth rate. In the study case, wave breaking and energy dissipation are most likely in autumn (October) between 44 and 52 km of altitude. Above this level, E p starts to increase again. These results are comparable with other Arctic, Antarctic and northern midlatitude lidar observations, as those of Alexander et al. (2011). With a vertical shear of stratospheric winds similar to our case, they find E p dissipation between 35 and 50 km and no dissipation in the autumn mesosphere. WRF seems to reliably reproduce the characteristics of summer and autumn atmospheric dynamics. Wave breaking mechanisms (by both wind filtering and saturation) are completely reproduced even if horizontal momentum fluxes, drag force and potential energy are underestimated with respect to the idealized experiment and lidar data.
The difference between real case and lidar data can be (at least) partly explained by the fact that WRF has no highly resolved information about (the numerous) thunderstorms occurring outside the inner domain. The amount of energy from such rain events can be partly or completely neglected. Energy is also underestimated because the spectrum of resolved small-scale waves is reduced. However, this is somewhat true also for lidar retrievals, with a band-pass filter centred at about 3.6 km of vertical wavelength, between 30 and 50 km, and 7.1 km above. We believe that filtering WRF wavelets to fit better with lidar spectral windows is not supposed to provide a valuable addition to the comparison. WRF energy is already underestimated and filtering would reduce this energy even more.
Further work is needed to analyse the sensibility of WRF to vertical and spatial resolution, and domain size, in order to find the right configuration that can ensure the best ratio between computational coast and realistic GW drag force and energy estimates. It is clear that a systematic comparison between model results and in situ data is a main way to achieve this issue. This would be the first step to use WRF as a fully complementary analysis tool, to ground-and space-based observations (limited to specific regions or certain latitude bands) to perform independent wavelet analysis and characterize convective GWs as a function of wave frequency. The implementation of a more accurate spectral parametrization, at different altitudes and latitudes, of HMF and drag force is a key point to improve significantly the reliability of global atmospheric circulation in weather forecasting and climate models (Butchart et al., 2010).