Impact of gravity waves on the middle atmosphere of Mars: a non-orographic gravity wave parameterization based on Global Climate modeling and MCS observations

The impact of gravity waves (GW) on diurnal tides and the global circulation in the middle/upper atmosphere of Mars is investigated using a General Circulation Model (GCM). We have implemented a stochastic parameterization of non-orographic GW into the Laboratoire de M\'et\'eorologie Dynamique (LMD) Mars GCM (LMD-MGCM) following an innovative approach. The source is assumed to be located above typical convective cells ($\sim$ 250 Pa) and the effect of GW on the circulation and predicted thermal structure above 1 Pa ($\sim$ 50 km) is analyzed. We focus on the comparison between model simulations and observations by the Mars Climate Sounder (MCS) on board Mars Reconnaissance Orbiter during Martian Year 29. MCS data provide the only systematic measurements of the Martian mesosphere up to 80 km to date. The primary effect of GW is to damp the thermal tides by reducing the diurnal oscillation of the meridional and zonal winds. The GW drag reaches magnitudes of the order of 1 m/s/sol above 10$^{-2}$ Pa in the northern hemisphere winter solstice and produces major changes in the zonal wind field (from tens to hundreds of m/s), while the impact on the temperature field is relatively moderate (10-20K). It suggests that GW induced alteration of the meridional flow is the main responsible for the simulated temperature variation. The results also show that with the GW scheme included, the maximum day-night temperature difference due to the diurnal tide is around 10K, and the peak of the tide is shifted toward lower altitudes, in better agreement with MCS observations.


Introduction
Gravity waves (GW) are small-scale atmospheric waves frequently detected in terrestrial planet atmospheres. They are an intrinsic feature of all stably stratified planetary atmosphere and play an important role in the large-scale circulation and variability of the middle/upper atmospheres of the Earth and Mars [Barnes, 1990;Joshi et al., 1995;Forget et al., 1999;Angelats i Coll et al., 2005;Fritts et al., 2006;Fritts and Alexander, 2003;Alexander et al., 2010;Medvedev et al., 2011]. Depending on the source, GW can be excited in the troposphere by a variety of mechanisms, which perturb the stratified atmospheric fluid and produce GW oscillations. The sources include flow over topography (orographic GW) or atmospheric convection, front systems and jet-streams (non-orographic GW). The restoring force is the buoyancy that results from the adiabatic displacements of air parcels characteristic of these disturbances [Fritts and Alexander, 2003]. GW propagation into the atmosphere provides a significant source of momentum and energy to the mean flow when they break or encounter critical values. Thus, the gravity wave drag (i.e. gravity wave momentum deposition to the mean flow) may cause wind acceleration or deceleration.
On the Earth, the observed reversal of the temperature gradient in the mesopause region, resulting in a winter polar warming around the mesopause level (∼ 80-90 km), is related to the effects of GW drag on the global circulation [Lindzen, 1981;Holton, 1982;Hauchecorne et al., 1987]. It is also well established that non-stationary (i.e. non-orographic) GW are a substantial driver of the quasi-biennial oscillation (QBO) in the equatorial regions of the Earth's atmosphere [Lindzen and Holton, 1968;Lott et al., 2012], complementing the forcing from the synoptic and planetary scale equatorial waves [Lott and Guez, 2013].
For Mars, several studies reported large density and temperature fluctuations (∼ 5 to 50 %) on small vertical scales height (<10 km vertically) from atmospheric entry profiles measured by Opportunity, Spirit, and Mars Pathfinder [Magalhães et al., 1999;Withers and Smith, 2006;Holstein-Rathlou et al., 2016], as well as small horizontal scales (typically 20-300 km) in various Mars accelerometer data sets at aerobraking altitudes [Keating et al., 1998;Withers, 2006;Creasey et al., 2006a;Fritts et al., 2006] and from mass spectrometer measurements and, Terada et al., 2017. Radio occultation temperature profiles by Mars Global Surveyor (MGS) obtained from the surface up to 35 km altitude also showed significant wave activity over the tropics, also partly attributed to zonally modulated thermal tides [Hinson et al., 1999;Creasey et al., 2006b]. Furthermore, "rocket dust storms" during which rapid and efficient vertical transport takes place by injecting dust particles at high altitudes in the Martian troposphere (30 to 50 km), may generate GW . The role of mesoscale GW is supposed to be crucial for local CO 2 condensation, responsible for the formation of mesospheric CO 2 clouds observed by Mars Express between 60 and 80 km altitude Montmessin et al., 2007].
MGS, Mars Odyssey (ODY) and Mars Reconnaissance Orbiter (MRO) aerobraking measurements were used to analyze density and temperature perturbations in the Martian thermosphere, where small-scale variability has been systematically observed whenever in situ data have been obtained (e.g. Keating et al. [2007]; Fritts et al. [2006]; Creasey et al. [2006a]; Zurek and Smrekar [2007]). MGS results (both radio profile and accelerometer data) leave open the question of the primary source mechanisms and source spectra, such as the distribution of phase speeds. Regarding the source, convection is often invoked to generate non-orographic GW (e.g. non stationary waves with non-zero phase speed) which propagate upwards Medvedev et al., 2011]. Fritts et al. [2006] and Creasey et al. [2006a] found that GW amplitudes vary significantly with time and season, being generally larger in Northern autumn and winter (i.e. second half of Martian year), and at middle to high latitudes, apparently reflecting mean source and filtering conditions. Their amplitudes also appear to vary with longitude and time and may provide clues to interactions with larger-scale motions [Fritts et al., 2006]. According to Creasey et al. [2006b], westwardpropagating waves may be encountering a critical level in spring and summer, but as the jet decreases towards autumn equinox, these waves may propagate to the thermosphere. More recently, NASA's Neutral Gas Ion Mass Spectrometer (NGIMS) instrument aboard the Mars Atmosphere and Volatile EvolutioN (MAVEN) satellite retrieved large (20-40 %) GW-induced CO 2 density perturbations in the Mars upper atmosphere between 180 and 220 km . Wave features were found with apparent wavelengths of ∼ 100 and 500 km in the Ar density around the exobase [Terada et al., 2017].
Many different parameterizations including those of Lindzen [1981]; Hines [1997] have been used to represent the impact of subgrid scale gravity wave processes in the terrestrial atmosphere. For Mars, the inclusion of orographic GW effects improves the performance of models at lower and higher altitudes [Barnes, 1990;Collins et al., 1997;Forget et al., 1999;Rafkin et al., 2001;Angelats i Coll et al., 2005]. Also, the thermal effect of non-orographic gravity waves has been proposed to explain some of the puzzling modelobservation discrepancy identified in the Martian atmosphere temperatures between 100 and 140 km [Medvedev and Yiǧit, 2012], although a major part of the discrepancy can be attributed to issues in the calculation of the non-LTE radiative cooling by CO 2 [Forget et al., 2009;Medvedev et al., 2015]. Medvedev et al. [2015] implemented a nonlinear spectral gravity wave parameterization [Yiǧit and Medvedev, 2010] in their GCM, extended up to 130 km, and their simulations showed that: 1. GW decelerate zonal winds at all seasons, 2. they produce jet reversals similar to those observed in the terrestrial mesosphere and lower thermosphere, 3. GW weaken the meridional wind and modify the zonal mean temperature by up to ± 15 K. However, modeling efforts for Mars suffer from lack of measurements to validate predicted wind fields and from no observational constraints on the GW forcing [Creasey et al., 2006a] Other studies performed high-resolution simulations (∼ 60 km grid size) with a general circulation model in order to resolve a significant portion of small-scale GWs and to capture the impact of GWs on the dynamics and energetics of Mars atmosphere without any parameterizations [Kuroda et al., , 2016[Kuroda et al., , 2019. Those authors showed that the GW activity varies greatly with season and geographical location, with stronger wave generation in the northern hemisphere winter in the mesosphere and smaller activity in polar regions of the troposphere throughout all seasons. Given the lack of global observations of GW the climatology of the small-scale disturbances obtained in Kuroda et al. [2019] can serve as a proxy for gravity waves that are largely not resolved by GCMs with conventional grid resolution.
In this paper we describe the results of the implementation of a non-orographic GW parameterization into the Mars Global Climate Model (MGCM) developed at the Laboratoire de Météorologie Dynamique (LMD), which already included an orographic GW scheme [Forget et al., 1999]. Although the LMD-MGCM was already able to reproduce temperature observations with an error of less than 10 K at most locations and times, we found that it was necessary to update the representation of several physical processes at work in the Martian atmosphere to further improve its accuracy. Among them: 1) the vertical distribution of the dust, characterized by detached layers (between 20 and 30 km) and by large day-night variation [McCleese et al., 2010;Heavens et al., 2011a,b,c;Madeleine et al., 2011;Navarro et al., 2014], 2) the radiative impact of clouds, which remains challenging to be well represented in the GCM , 3) the effect of GW (orographic and/or non-orographic), 4) the phasing of the thermal tide wave in the vertical compared to data from the Mars Climate Sounder (MCS) on board MRO . MCS-MGCM biases probably resulted from the incorrect representation of the three previous processes mentioned above  and several improvements are currently being implemented.
Here we took advantage of the previous development work carried out by our team for the LMD Earth [Lott et al., 2012] and Venus GCM  to also implement the same non-orographic GW scheme in the MGCM. According to Lott et al. [2012], the stochastic parameterization represents the unpredictable aspects of the sub-grid dynamics, considering that convection and fronts can generate waves throughout the full range of phase speeds, wave frequency, vertical and horizontal scales.In a context where the sources are not known precisely, we consider this is a reasonable choice. One of the main goals of this study is to understand the impact of non-orographic GW on the global circulation and the thermal structure of the Martian middle atmosphere (50-100 km altitude). What is the magnitude of GW-induced drag? Where do GW break/saturate and deposit the maximum momentum? Can GW explain the remaining discrepancy between the MCS observations and the GCM simulations? If so, what is their impact on the predicted winds?
Section 2 describes the non-orographic GW scheme and the main tunable parameters adopted in this work. The LMD-MGCM model and MCS dataset used here are described in Sections 3 and 4, respectively. The impact of our non-orographic GW scheme on the simulated wind and temperature is described in Section 5. The method implemented here to identify the set of "best-fit" GW parameters is explained in Section 6, together with the comparison with MCS results. A number of sensitivity tests is discussed in Section 7 and concluding remarks are given in the Section 8.

Non-orographic GW parameterization: formalism
The scheme implemented here is based on a stochastic approach, as fully described in Lott et al. [2012] and Lott and Guez [2013]. In such a scheme a finite number of waves (say M = 8), with characteristics chosen randomly, are launched upwards at each physical time-step (δt = 15 min) and at each horizontal grid point.The question of the location of the non-orographic GW sources is a difficult one and given the lack of information on such location, launching the gravity waves at each grid point is probably the simplest assumption. The number of waves is given as M = NK × NO × NP, with NK=2 values of GW horizontal wave-numbers, NO=2 absolute values of phase speed, and NP= 2 directions (westward and eastward) of phase speed. This approach allows the model to treat a large number of waves at a given time t by adding the effect of those M waves to that of the waves launched at previous steps, to compute the tendencies. We need to parameterize GWs whose life cycle (i.e. from its generation to wave break) is contained in a characteristic time interval ∆t that has to significantly exceed the GCM time step δt. On the Earth, GW theory indicates that atmospheric disturbances induced by convection have life cycles with duration ∆t around 1 day (∆t = 24h) [Lott and Guez, 2013]. This typical time scale is likely to be relevant for Mars too: the choice for this timescale on Earth is made considering a mix of inertio-gravity waves (influenced by the planetary rotation, which is similar on both planets) and convectivelygenerated gravity waves, which have shorter frequency than the inertio-gravity waves, but which can propagate for several hours before dissipation or breaking occurs. The spectrum is discretized in 770 stochastic harmonics (≈ M x ∆t/δt) which contribute to the wave field each day and at a given horizontal grid point. At each time t the vertical velocity field w of an upward propagating GW can be represented by this sum: where C n are normalization coefficients such that ∞ n=1 C 2 n = 1. It is assumed that each of the w n can be treated independently one from the others, and each C 2 n can be viewed as the probability that wave field is given by the GW w n (see also the expression for C n in Equation  6). This formalism is then applied to a very simple multi-wave parameterization, in which w n represents a monochromatic wave as follows where the wavenumbers k n , l n and frequency ω n are chosen randomly. In Equation 2, H∼ 11 km is a middle atmosphere characteristic vertical scale for Mars and z is the logpressure altitude z = H ln(P r /P), with P r a reference pressure (P r = 250 Pa), taken here at the height of the source (i.e. above typical convective cells z ∼ 8 km). To evaluate the amplitude ofŵ n we randomly impose it at a given launching altitude z 0 , and then iterate from one model level z 1 to the next z 2 by a Wentzel-Kramers-Brillouin (WKB) approximation (see Equation 4 of Lott et al. [2012] for details). Using that expression plus the polarization relation between the amplitudes of large scale horizontal windû and vertical windŵ (not shown here), we can deduce the Eliassen-Palm (EP) flux (vertical momentum flux of waves), with k, l the horizontal wavenumber and ω the frequency of the vertical velocity field. The latter is included in the vertical wavenumber m = N | ì k | Ω taken as the WKB non-rotating approximation in the limit H → ∞, with Ω = ω − ì k ì u and N the Brunt-Vaisala frequency (see Lott et al. [2012] for details). In Equation (3) ρ r is the density at the reference pressure level P r (ρ r ∼ 0.007 kg m −3 ).
In our scheme we randomly impose bothŵ n and the EP-flux at a given launching altitude z 0 (see Section 2.1). We also have explicitly introduced a constant vertical viscosity µ (see Eq. 4 in Lott et al. [2012]) which controls the GW drag vertical distribution near the model top. To move from one model level to the next model level above, we essentially conserve the EP-flux, but allow a small diffusivity, ν = µ/ρ 0 , which can be included by replacing Ω by Ω + iνm 2 . This small diffusivity is here to guarantee that the waves are ultimately dissipated over the few last model levels, if they have not been before (hence the division by the density ρ 0 ). In addition, this new EP-flux amplitude is limited to that produced by a saturated monochromatic waveŵ s following [Lindzen, 1981]: or eitherŵ = 0 when Ω changes sign, to treat critical levels. In Equation 4 S c is a tunable parameter and k * a characteristic horizontal wavelength corresponding to the longest wave being parameterized (see more details in Sec. 2.1) Finally, we calculate the tendencies ρ −1 δ z ì F z n (n = 1, M) produced by the GW drag on the winds, computing the tendencies due to the M generated waves. Since we assumed that w n are independent realizations, the mean tendency they produce is the average of these M tendencies. Thus, we first redistribute the averaged tendency over the longer time scale ∆t by re-scaling it by δt/∆t and second, we use the auto-regressive (AR-1) relation described in [Lott et al., 2012] as follows: This indicates that, at each time step, we promote M new waves by giving them the largest probability to represent the GW field, and degrade the probabilities of all the others by the multiplicative factor (∆t − δt)/∆t. As explained in Lott et al. [2012], by expressing the cumulative sum underneath the AR-1 relation in Equation 5, we recover the formalism for infinite superposition of stochastic waves by taking: where p is the nearest integer that rounds (n − 1)/M (i.g. toward lower values).

GW parameters setup
In this section we describe the main tunable parameters used in the non-orographic GW scheme implemented here. The characteristics of every wave launched in the GCM are selected randomly with a prescribed box-shaped probability distribution, whose boundaries are key model parameters. These are chosen on the basis of observational constraints (whenever available) and theoretical considerations, and tuned using MCS data comparison. A total of 26 runs for the full Martian year (MY) 29 have been performed to cover a large range of possible combinations of the GW characteristics, based on those constraints. One of the advantages of the GW scheme used here is that, contrarily to other GW parameterizations, each parameter has a physical meaning, as described below.
Source height and duration First, we assumed that the non-orographic GW source is placed above typical convective cells (i.e. around 8 km, depending on the topography). Mars Express Radio Occultation (RO) profiles provided good coverage at latitude and local times where planetary boundary layer (PBL) convection is occurring, giving an accurate determination of the depth and the spatial variation of the convective boundary layer (CBL) at fixed local time (∼ 17 h) [Hinson et al., 2008]. Those authors found that the CBL extends to an height of 3-10 km above the surface at the season of the measurements (mid-spring in the northern hemisphere). Spiga et al. [2010] compared the RO temperature profiles to large-eddy simulations performed with the LMD Martian mesoscale model and found intense CBL dynamics within the measured depths (up to 9 km). In our scheme the parameter controlling the wave launching level z 0 is σ = P/P s = 0.4, that corresponds to a region centered at about 250 Pa (∼ 8 km), covering pressures from 160 Pa to 348 Pa, depending on the surface pressure P s and varying with season. The source is chosen to be uniform, without latitudinal variation, and aside from the atmospheric profiles through which the GWs' propagation is modeled, nothing else in the scheme varies with location or season. Considering that non-orographic GW are expected to be generated by multiple sources (e.g. PBL convection, jet acceleration, dusty convection etc) and that represents a complex interplay of timescale for GCMs, we assume here that the source is turned on all day. In future developments, we foresee to implement in our scheme a characteristic time interval more representative of the life cycle of GW produced by PBL convection on Mars, but this is beyond the scope of this paper. EP-flux amplitude F z from Equation 3 gives the vertical rate of transfer of wave horizontal momentum per unit of area. This value has never been measured in Mars' atmosphere and it represents an important degree of freedom in the parameterization of gravity waves. Thus, in our scheme we impose the maximum value of the probability distribution F z max at the launching altitude z 0 (see previous point), for every set of MY29 runs. In order to define F 0 max we have explored typical values used in the literature to evaluate the order of magnitude of the EP-flux within the realm of what is realistic. For instance, using aerobraking data Fritts et al. [2006] found an estimation for momentum flux about 2000 m 2 s −2 at altitudes 100-120 km where the density is typically from 10 −7 to 10 −9 kg m −3 , and this yields values for EP-flux of the order 2·10 −6 to 2·10 −4 kg m −1 s −2 . Calculations by a one-dimensional full-wave model performed by Parish et al. [2009] gave mean momentum flux (e.g.ρu w ) ≈ 7·10 −7 kg m −1 s −2 below 100 km, at 82 • N latitude during winter solstice, for a GW packet propagating with horizontal wavelength between 38 km and 150 km (see their Figure  6). Other authors [Medvedev et al., 2011 also implemented a non-orographic gravity wave scheme [Yiǧit and Medvedev, 2010;Yiǧit et al., 2015] in their GCM using an approach different from ours. They adopted an analytical form of the GW spectrum, in which the vertical propagation of horizontal momentum fluxes is given by , withū 0 the mean wind at the source level, c j the phase speed of the harmonic j, and c w the half-width at half maximum of Gaussian distribution, fixed at c w = 35 m s −1 . In their experiments the spectrum is discretized with 28 harmonics for two values ofū 0 = 0 and 20 m s −1 . Although they employed a geographically uniform flux per unit mass u w max = 2.5·10 −3 m 2 s −2 , the wave source varies in time and space because of the modulation by the simulated local windū 0 in the lower atmosphere. In their scheme the direction of propagation of GW harmonics coincides with that of local wind at z 0 for c j > 0 (positive momentum fluxes) and it is against it for c j < 0 (see Medvedev et al. [2011] for more details.) Assuming typical mean densities from 10 −3 to 10 −2 kg m −3 around 250 Pa (i.e. the GW source level of our simulations) this gives a momentum flux F 0 max above typical convective cells of the order of 10 −6 to 10 −5 kg m −1 s −2 , which is in the range of values explored in our study. It should be stressed here that the stochastic approach implemented in our scheme has the advantage of allowing to treat a wide diversity of emitted gravity waves, thereby a wide diversity of momentum fluxes. Our only setting is the maximum EP-flux amplitude at the launching altitude, therefore we tested several values of F 0 max by performing sensitivity tests for values of F 0 max ranging from 10 −10 kg m −1 s −2 to 10 −5 kg m −1 s −2 (see also Section 7). Then the amplitude of the EP-flux for each GW is chosen randomly between 0 and F 0 max . We set 10 −10 kg m −1 as a lower limit for our tests because GW impact is negligible in the LMD-MGCM for momentum flux smaller than that value. Horizontal wavenumber In our scheme the horizontal wave number amplitude is defined as in Lott et al. [2012] k * < |k | < k s . The minimum value is k * = 1/ √ ∆x∆y, where ∆x and ∆y are comparable with the GCM horizontal grid (δx and δy ≈ 600 km in this study). This value corresponds to the longest waves that one parameterizes with the current GCM horizontal resolution. The maximum (saturated) value is k s < N/u, N being the Brunt-Vaisala frequency associated to the mean flow and u 0 the mean zonal wind at the launching altitude. The corresponding min/max values for the GW horizontal wavelength (λ h = 2π/k ) are between 10 km and 300 km. Those values are within the observed range of GW wavelengths [Magalhães et al., 1999;Hinson et al., 1999;Fritts et al., 2006]. Phase speed Another key parameter is the amplitude of absolute phase speed |c| = |ω/k |.
As for the other tunable parameters, we impose the minimum c min and maximum c max values of the probability distribution at the beginning of the runs, and the model chooses randomly |c| between c min and c max . Here c min is set to 1 m/s (i.e. for nonstationary GW) and c max is of the order of the zonal wind speed at the launching altitude. The range of c max is consistent with previous values used in the literature . Both eastward (c > 0) and westward (c < 0) moving GWs are considered. Three values of maximum probable c max have been tested to evaluate the sensitivity of winds and temperature to this parameter (i.e. 10 m/s, 30 m/s and 60 m/s). See Section 7 for details. Saturation parameter S c is a tunable parameter in our scheme, on the right hand side of Equation 4, which controls the breaking of the GW by limiting the amplitude w s . In the baseline simulation we set S c = 1.

LMD-MGCM description
The LMD-MGCM is a finite-difference model based on the discretization of the horizontal domain fields on a latitude-longitude grid [Forget et al., 1999]. The horizontal resolution used in this work is 64 longitude x 48 latitudes (3.75 • x 5.62 • ). Vertical levels are hybrid coordinates, with 32 levels up to 0.05 Pa (about 100 km altitude), above the top of MCS profiles. This is a standard vertical resolution, as used in all previous recent papers describing LMD-GCM and MCD results  and . The GW parameterization is precisely designed to cope with the coarse vertical resolution of the GCM: the propagation of the GW is not resolved by the GCM, it is a subgrid-scale phenomenon whose impact on the large-scale flow (heat and momentum transfer when they break) is computed by the parameterization, then passed to the dynamical core of the GCM as a tendency for temperature and wind. For each physical timestep of 15 Martian minutes (1 minute is about 1/1440th of a Martian sol), 10 dynamical time steps occur. The model includes the CO 2 cycle, with condensation and sublimation of CO 2 and the global change of the total atmospheric mass [Forget et al., 1998]. The cycles of dust and water and their interactions are also represented. Dust is modeled with a two moments scheme, that transports both the mass mixing ratio and the number of dust particles, with a distribution for particle size assumed to be log-normal with a fixed effective variance [Madeleine et al., 2011]. Lifting is global and constant with a rescaling applied to the total column quantity at each physical timestep to match the observed column of dust during Martian Year 29 [Montabone et al., 2015]. The water cycle includes a microphysical scheme for the sublimation and condensation of water ice clouds on the dust particles and interactions between the ice at the surface and the atmosphere . Dust and water ice in the atmosphere are 4D variables, transported and radiatively active [Madeleine et al., 2011[Madeleine et al., , 2012. Another key improvement was the parameterization of convection and near surface turbulence using a thermal plume model, which is coupled to surface layer parameterizations taking into account stability and turbulent gustiness to calculate surface-atmosphere fluxes [Colaïtis et al., 2013].
The LMD-MGCM also includes a parametrization of the orographic gravity waves [Forget et al., 1999;Angelats i Coll et al., 2005] based on the gravity wave drag scheme of Lott and Miller [1997] and Miller et al. [1989]. At the time the scheme parameters were chosen conservatively and the effect of the orographic wave drag was probably underestimated. In this work, we add a more general scheme which includes non-orographic gravity waves, but it is possible that this also accounts for the previously underestimated effect of orographic gravity waves. At this stage it is very difficult to separate the contribution of orographic sources to the wind drag from non-orographic ones, but this will be certainly the next step in a future theoretical work. MCS [McCleese et al., 2007] is a limb-viewing infrared radiometer aboard the sunsynchronous MRO spacecraft. Its nominal local times of observations at most latitudes (except when the orbit crosses high latitudes) are fixed around 03:00 and 15:00 h Local Time (LT). Kleinböhl et al. [2009] obtained profiles of temperature (using the CO 2 15 µm absorption band) as well as profiles of dust and water ice extinction opacity (at wavelengths centered respectively around 21.6 µm and 11.9 µm), nominally from the surface to about 0.1 Pa with 5 km vertical resolution. MCS data [McCleese et al., 2010] currently provide the only systematic measurements of temperature in Mars' mesosphere up to about 80 km. For several technical reasons explained in the above-cited papers, though, a large number of profiles stop well before reaching the surface. Although the nominal vertical resolution of the instrument is 5 km (corresponding to the separation of the peaks of the weighting functions), the profiles are over-sampled using information from more than one weighting function. Both temperature and aerosol extinction opacity profiles are standard MCS products. Most recently, Kleinböhl et al. [2017] developed a new scheme for the retrievals of these quantities, based on a two-dimensional rather than one-dimensional radiative transfer. Instead of simply assuming spherical symmetry, the 1-D retrieved fields are interpolated along the orbit to provide correction factors. Further retrievals are then performed using the interpolated fields and a 2-D radiative transfer scheme, in order to correct for gradients along the line-of-sight. This new retrieval methodology (producing retrieval version 5.2) is particularly useful to mitigate dayside-nightside differences in retrievals at high latitudes, where horizontal gradients may be strong.

MCS Data set description
In this paper, we use these latest retrievals (version 5.2). Temperature uncertainties are typically around 0.5 K and only increase at low altitudes (below about 5 km altitude), where the atmosphere starts to become opaque, and at altitudes above about 60 km, where the instrument signal-to-noise ratio starts to decrease. Since we are interested in comparing gridded GCM temperature fields with data, we have carried out binning of MCS observed temperatures in 4-D bins (latitude, longitude, solar longitude, and local time). We separate observations into dayside (06:00 h < LT ≤ 18:00 h) and nightside (00:00 h < LT ≤ 06:00 h and 18:00 h < LT ≤ 24:00 h) bins. As mentioned above, most MCS observations at latitudes equatorward of |80 • | are provided around 03:00 h in the nightside and 15:00 h in the dayside. Exceptions are cross-track observations that span additional local times, but this kind of observation only started after September 13, 2010 (L S = 146 • , MY 30) [Kleinböhl et al., 2013], therefore it does not affect the results of this paper, in which we only use MY29 observations. This particular year was chosen because at the time of the preparation of the manuscript it was the most complete year after the global dust storm in MY28. A comparison with less dusty years (but with worst coverage) such as MY30 or MY31 is left for a future study. Furthermore, we limit our comparison between model results and data to latitudes equatorward of |80 • |. By doing this, we can safely use GCM results at local times 03:00 h in the nightside and 15:00 h in the dayside, and avoid comparison at non-homogeneous local times. The width of the latitude, longitude, and solar longitude bins of MCS observations is, respectively, 3.75 • , 5.625 • , and 5 • , in order to match MGCM horizontal resolution. The binned values are provided on the same pressure grid as the one used by the MCS team for their standard products. It is worth noting that MRO entered a long period of safe mode in late 2009, therefore there are no MCS observations in MY 29 after L s = 328 • .

Impact of non-orographic GW drag on LMD-MGCM zonal winds
Examples of simulated zonal mean winds before and after the implementation of the non-orographic GW parameterization in the LMD-MGCM are plotted in Figures 1-4 at four solar longitudes: L s =0 • , L s =90 • , L s =210 • and L s =270 • . The runs described here correspond to the MGCM with the GW scheme on, using the reference parameters listed in Table  1. The choice of this set of parameters will be discussed in section 6.
Without the non-orographic GW scheme activated the simulated wind structure is comparable to the results described in previous theoretical works [e.g. Forget et al., 1999;Haberle et al., 1999;Hartogh et al., 2005]. It shows a strong (retrograde) easterly jet at northern hemisphere spring equinox (L s =0 • ), in the middle atmosphere above ∼ 0.5 Pa (about 60 km) at equatorial regions, and strong (prograde) westerly wind in both hemispheres at high latitudes with velocities increasing with height (see panel a in Figure 1). At the northern winter solstice the circulation of Mars atmosphere is dominated by a quasi-global Hadley cell extending from 60 • S almost to the north pole and model simulations predict a strong prograde wind jet corresponding to the pole-to-pole heating gradient (panel a in Figure 4). When the atmospheric dust loading has one of its peaks (L s =210 • ), the ascending branch of the Hadley cell reaches higher latitudes if the thermal forcing is very strong, for instance when the dustiness increases, and the prograde winds are also stronger than at other seasons (panel a in Figure 3). At northern summer solstice instead (panel a in Figure 2), MGCM simulations show a weaker thermal forcing than at northern winter solstice, because of the smaller amount of dust in the atmosphere and the reduced solar flux near aphelion [Forget et al., 1999]. Consequently, the Hadley circulation is much less intense and the simulated southern winter polar warming also weaker than at northern winter solstice.
When the non-orographic GW scheme is on, the retrograde wind jet would prevent most of the westward (c < 0) propagating GW from reaching the mesosphere, leaving primarily eastward waves with positive momentum flux to propagate, with the reverse being true for westerly jets at mid-high latitudes. The complete wave momentum deposition occurs in those layers where upward propagating non-orographic GW energy with phase speed c meets a zonal wind vector opposite to its direction. In other words, when c is close to the mean zonal wind (c →ū) or if its amplitude is large enough to create shear instability (c −ū →ū 0), the saturation of a gravity wave occurs and the wave momentum flux is transferred to the mean flow [Lindzen, 1981;Terada et al., 2017]. For non-orographic GW with a given phase speed c ū, the intrinsic speed |ū − c| becomes larger in amplitude in presence of wind jets, inducing a stronger GW drag (i.e. deceleration if c <ū or acceleration c >ū). Therefore GW drag is especially strong in the upper part of the zonal winter jets, where the contrast between the mean flow velocity and the given phase speed (toward which the speed of the flow is decelerated/accelerated) is high. The induced drag due to the non-orographic GW momentum deposition to the mean state changes the large-scale horizontal winds, hence the large-scale circulation and vertical winds. Consequently, those waves drive large-scale heating/cooling by adiabatic downward/upward large-scale vertical motion associated with the altered circulation.
The GW parameterization scheme implemented here allows for quantifying the "GW drag" (e.g. the zonal and time averaged gravity wave momentum deposition to the zonal flow) in order to gain further insight into the nature of changes in the modeled mean fields (see panels b in Figures 1-2). In all simulations the GW drag is significantly large (> |0.001| m/sol/s) above ∼ 1 Pa (about 50 km) and it predominantly reduces the westerly jets at high latitudes (negative drag) and weakens the equatorial easterly jet by accelerating the zonal wind (positive drag). In most seasons the negative GW drag on the retrograde zonal wind is about two orders of magnitude larger than the positive one for the same period. This is due to the damping of diurnal tides by the GW (see Section 6.2), and thus of the mean zonal wind: the diurnal tides tend to drag the zonal wind toward the phase velocity of the tides (i.e. the motion of the sub-solar point in a westward direction) which is about -240 m/s, and the retrograde jet is mainly forced by the interaction between thermal tides and wave mean flow [Forbes et al., 2002]. During the dust season (L s =210 • ) westward zonal wind jets are stronger than at other seasons above 1 Pa, therefore the impact of GW on the general circulation is expected to be larger. Our results in Figure 3 show instead that the GW-induced deceleration during this season is similar to that in Figure 4 (L s =270 • ) at same pressure, reaching approximately 0.55 m/s/sol between 10 −1 Pa and 10 −2 Pa. However, those levels are close to the upper limit of our runs, and further studies with a vertically extended version of the LMD-MGCM should be performed in order to better quantify the GW drag at those altitudes.
A tentative comparison with previous works is give here. Unfortunately, there are very few published studies on non-orographic GW parameterization implemented into Mars GCMs, and they are mostly focused on the theoretical impact of GW in the thermosphere [Medvedev et al., 2011;Medvedev and Yiǧit, 2012]. Those authors also claimed that the main influence of GW in the winter hemisphere is near the edge of westerly jet, and that the net effect on the mean zonal wind is to decelerate and even reverse it increasingly with height (between 100 and 130 km). The magnitudes of the drag they found vary from tens at 0.1 Pa (∼ 80 km) to hundreds of m/s/sol above, depending on the shape of simulated jets, and the assumed wave sources. In our simulations GW drag reaches maximum magnitudes of the order of 1 m/s/sol around 10 −2 Pa in the northern hemisphere (NH) winter solstice, about a factor 100 smaller than values estimated by Medvedev et al. [2011] at the same season (see their Figure 1). This discrepancy may be related to the difference between the maximum amplitude of the EP-Flux at the launching altitude used in our best-fit simulations (F 0 max = 7 · 10 −7 kg m −1 s −2 ) and the EP-flux values used in Medvedev et al. [2011], which are about 2 order of magnitude larger (10 −4 -10 −5 kg m −1 s −2 ). Our maximum value is however consistent   Parish et al. [2009] and with the lower limit estimated by Fritts et al. [2006]). Nevertheless, a quantitative comparison of GW-induced wind drag should be taken with caution, because it would require first a detailed inter-comparative study between the two different GW parameterizations, which is not straightforward and should be addressed in a future dedicated paper. Here we have investigated the possible impact of the vertical resolution on the GW drag by increasing the vertical resolution from 6-7 km (32 levels) to 2-3 km (54 levels). The results (see Figure in Supplement Material) show that GW drag values are of the same order of magnitude with higher resolution than with our standard GCM model resolution, confirming that our GW parameterization does not depend on the vertical resolution."

Comparison with wind observations
As shown in the previous section, the non-orographic GW parameterization implemented in this work produces major changes on the simulated zonal winds. Overall, the GW slow down zonal winds everywhere, and in the tropics they become nearly zero (panels c  in Figures 1-2). Unfortunately, there are no direct wind observations to compare with, except for rare and sparse observations from ground based telescopes [Sonnabend et al., 2012;Lopez-Valverde et al., 2016]. The wind fields derived from MCS temperature fields using geostrophic balance are based on several assumptions (see McCleese et al. [2010]) and they are certainly not the gold standard to validate model simulations. They are estimates of zonal gradient wind calculated via thermal wind, but we can expect the atmosphere of Mars to be more complicated than that. Here instead we compare our results with the retrieved wind velocity in the mesosphere of Mars during one observing campaign occurred in MY29 zonal wind ± 30 m/s, where a MGCM with zero wind in that region is consistent. The same is true for the campaign A occurred in MY30 (not shown here), where low westerlies winds around 20-30 m/s at southern latitudes are also consistent with values simulated in this work. In northern latitudes the fit is not very good (both with or without GW) but there the data is also all over the place, with 30 m/s westward in campaign A and 180 m/s eastward in campaign B, at a similar season, and measured uncertainties are also larger.

Impact of GW mean flow forcing parameters
The momentum deposited by GW on the mean state not only may reverse easterly solsticial equatorial zonal jets in the middle atmosphere (∼ 5-0.05 Pa) and weaken the west-   erly zonal jets at mid-high altitudes, as described in section 5, but it also drives changes in the meridional flow. Figures 6, 7, 8 show the alteration of zonal wind, meridional circulation and temperature due to parameterized GWs, and the GW drag on zonal wind for the northern hemisphere (NH) Spring (L s =0,30 • ) , NH summer (L s =90 • ,120 • ) and NH Winter (L s =270 • ,300 • ), respectively. Note that the contour maps in panels c of Figures 6-8 represent monthly averaged GW drag, while Figures 1, 2 and 4, show daily averaged at the corresponding solar longitude L s = 0 • , L s = 90 • and L s = 270 • , respectively. Also note that the contour scales are different from Figures 1-4. Focusing on Figure 6, the equator-to-pole flow increases by ten m/s at about 0.5 Pa (60 km approximately) with respect to the runs without GW (see panel f), thus increasing convergence at the pole, hence downward motions. As a consequence, the polar warming increases up to 30 K in both hemispheres, and the adiabatic cooling increase up to 10 K in the layers above, at mid-latitudes. At the equator, local heating above 10 −1 Pa (around 80 km) results from the reduction of ascending branch of the Hadley cell at these altitudes. For the period L s =90 • ,120 • (Figure 7), the dynamical mechanism is analogous, but the impact is smaller: the southern polar warming is increased by more than 15 K, caused by the GW induced acceleration of the north-to-south pole meridional flow above 1 Pa (see panel f). Similarly, in the NH winter ( Figure 8) the effect of GW is to accelerate the upper branch of the solsticial Hadley cell up to 30 m/s, and adiabatic heating takes place due to the intensification of poleward circulation cells, thus increasing the northern polar warming in the middle atmosphere of Mars by about 15 K on average. GW drag reaches magnitudes of the order of 1 m/s/sol above 10 −2 Pa in the NH winter solstice (panel b in Figure 2), and produces a major change in the zonal wind field (∼ 100 m/s), while the impact on the temperature field is relatively moderate (∼ 10-20 K). Those considerations indicate that GW induced alteration of the meridional flow is responsible for the simulated temperature variation. The impact on the temperature field will be also discussed in the next session.

Results: comparison with MCS data
Since the characteristics of GW spectrum are not well known on Mars, the strategy adopted in this work was to identify a set of reference tunable GW parameters that reduced model data biases by the greatest amount by comparing the LMD-MGCM and MCS thermal structure, specifically using diurnal tides as diagnostics. We run 26 simulations of a full Martian Year (MY 29) with non-orographic GW parameterization included. The subset of "bestfit" GW characteristics listed in Table 1 was selected with the help of sensitivity tests (see Section 7). Both MCS data and LMD-MGCM simulations were binned in boxes of 3.75 • latitude x 5.625 • longitude.
In this section we will focus only on the LMD-MGCM simulations performed with the non-orographic GW scheme activated using the wave characteristics as in

Impact of non-orographic GW drag on LMD-MGCM thermal structure
Examples of the impact of our GW scheme on the thermal structure of the Martian mesosphere predicted by the LMD-MGCM (up to about 100 km) are given in Figures 9-18. Figures 9-13 compare zonal mean day-night temperature averages T mean from MCS data (top panels) with model results from the LMD-MGCM with and without the non-orographic GW parameterization (middle panels). Differences between the model runs with and without the parameterization are also given. T mean stands for (T am + T pm )/2, where T am and T pm are the temperature at 03:00 h and 15:00 h local time, respectively. As discussed in Section 4, latitudes larger than 80 • (North and South) are excluded from this study to avoid comparison at non-homogeneous local time. In Figures 17 and 18 we plot examples of nightime (T am ) and daytime (T pm ) zonal mean temperatures at pressure z= 1 Pa, for latitude ranges 20 • S-20 • N and 60 • N-80 • N, respectively, during MY29 and for all seasons in which MCS data are available (L s = 0-330 • ).
It is assumed that the state of the model before including the new GW parameterization was consistent but incomplete to represent the reality, and we add the new GW parameterization hoping to get a more realistic emulation of the Martian atmosphere by our GCM. However, there is a chance that the new parameterization is instead compensating for errors in other routines. This weakness applies to any GCM study and we have to make the underlying assumption that there is no unphysical error in the existing physics routines implemented in our model. Notably, the parameterization by Madeleine et al. [2014] and Navarro et al. [2014] indeed improved the comparison between the LMD-MGCM and the observations for good physical reasons (e.g. coupled dynamical-radiative processes for dust, and radiative effect of clouds), but we have no choice other than to conduct our GCM work under the above-mentioned underlying assumption.

Equinoxes
The general circulation at the Mars equinoxes is dominated by two prograde midlatitude jets corresponding to the equator-to-pole heating gradients at low altitude with a single Hadley cell in each hemisphere. Above 10 Pa (∼ 40 km), the LMD-MGCM predicts a dynamically driven temperature inversion around 1 Pa (∼ 60-70 km) above both poles [Forget et al., 1999]. When the GW scheme is off, we found that during the Northern Hemisphere (NH) Spring Equinox (see Figure 9) the MGCM-MCS difference exceeds 20 K; in particular the model is colder at middle to high latitudes below 1 Pa, and warmer in the region between 1 and 0.1 Pa (about 50-75 km altitudes), mainly around 60 • latitude (North/South). In the MGCM version with the non-orographic GW scheme on, the region between 50 • and 80 • latitude (North/South) is up to 22 K warmer below 1 Pa and cooler above between 40 • and 70 • latitude (North/South), thus reducing the differences with MCS data. The bottom panel in Figure 9 illustrates the net effect of our non-orographic GW parameterization on zonal mean average temperature. Temperature in the tropical region (40 • S-40 • N) is also reduced by about 8 K around 1 Pa, in better agreement with the data. Similar conclusions hold when comparing results for NH Autumn (Figure 10). However, in this case the warmer region around 1 Pa around the equator is only partially reduced when the non-orographic GW scheme is activated. In Figures 9 and 10, the impact of GW on the temperature field is the one described for orographic (i.e. low phase speed) GW on Mars by early modeling studies like Barnes [1990]; Collins et al. [1997]; Forget et al. [1999]: the adiabatic cooling at low latitudes and the polar warming is enhanced and shifted to lower altitudes (10 to 100 Pa level) as a results of the GW friction acting on the zonal wind. Lowering the zonal wind reduces the Coriolis forces that limit the poleward meridional winds. In particular, this enhances the mass convergence at high latitudes and strengthens the polar warming.

Solstices
Concerning the solstice, the impact of non-orographic GW in MGCM simulations is also significant, during both the NH summer solstice (L s =90 • -105 • ) and NH winter (L s =270 • -285 • ), as shown in Figure 11 and 12, respectively. As discussed in Forget et al. [1999], around Northern winter solstice the strong pole-to-pole diabatic forcing creates a quasi-global Hadley cell which extends to 0.05 Pa (∼ 80 km). In such a cell the Coriolis force contributes to accelerate the poleward meridional motion on Mars, thus inducing a mass convergence and strong warming of the middle polar atmosphere down to about 5 Pa (∼ 25 km), as also observed [Jakosky and Martin, 1987;Theodore et al., 1993]. Forget et al. [1999] suggested that thermal inversions can generally be expected around 1 Pa (60-70 km) above the winter polar regions near solstice, and above both poles near equinox. However, we note here that the effect of GW is not an enhancement of high latitude (60 • -80 • ) warming (see also Figures 14  and 13), as during the equinoxes. On the contrary, around Northern winter solstice the GW drag tends to reduce the high latitude warming by more than 15 K (e.g. Figure 12), while heating the atmosphere above 1 Pa (∼ 60 km) at low/middle latitudes. This appears different than in most published results on orographic and non-orographic GW, and the opposite of what was simulated at the equinoxes (compare bottom panels in Figures 9 and 12). In reality, as discussed in Section 5.2, the polar warming (poleward of 80 • N, not shown in Figures 12) is increased as expected (see Figure 8), and it is only between 30 • N and 80 • N that the warming is reduced. This results from the poleward shift of the warming subsidence, but also from the direct effect of the GW drag on the meridional wind and on the damping of thermal tides (see section 6.2), which partly controls the high latitude warming around northern winter solstice [Wilson and Hamilton, 1996]. Data-model comparison improves mostly in NH winter, with biases at pressure range 10-1 Pa reduced up to 20 K, especially between 40 • N and 80 • N latitude (see Figure 12). There are still differences between MGCM and MCS results: at similar pressure levels and around the equator simulations are warmer than MCS data by 5 K, thus increasing model biases. Possible causes for those remaining discrepancies are related to: i) the uncertainties on the dust distribution in the model during the dust season, ii) the radiative effects of water-ice clouds on the thermal structure during the NH summer, and/or iii) the vertical resolution in the water ice cloud structure, during the northern summer season. At the aphelion the uncertainties related to the cloud modeling are larger than the effect of non-orographic GW on the large-scale circulation, and it is still challenging to reach a good accuracy to represent it well with GCMs.

Thermal structure seasonal variations
To illustrate the seasonal evolution of simulated temperature after the implementation of the GW scheme into the LMD-MGCM using the baseline GW parameters as in Table 1, zonal day-night averages (T am + T pm )/2 are plotted in Figures 13-16 as function of solar longitude L s for a selection of latitudinal bands: 20 • S-20 • N, 40 • N-60 • N, 60 • N-80 • N and 60 • S-80 • S. Overall, those figures indicate that there is a small reduction in the magnitude of the biases at most seasons and latitudes, with some notable exceptions: (i) the midaltitude negative bias at high southern latitudes in the first half of the year (northern spring and summer) is greatly reduced up to 20 K (Figure 13), (ii) the mid-altitude positive bias at high northern latitudes in early winter (L s = 240-270) also decreases by 15-20 K (Figure 14). Note that this period corresponds to the bulk of the dust storm season, yet the GW scheme improves the predicted thermal structure at those latitudes (iii) Positive biases in northern mid-latitudes (40 • -60 • N) are generally reduced up to 10 K at all altitudes and times of the year, except at very high altitudes around local winter solstice (Figure 15) (iv) Large negative biases at low altitudes at high latitudes (in northern fall) and northern mid-latitudes (in northern fall and winter) are not improved by the non-orographic GW parameterization and may be due to incorrect dust distributions and/or interactions between dust and water ice clouds during the dustiest time of the year. Remaining biases in the dust storm season (L s = 180 • -330 • ) between GCM simulations with the GW scheme activated and MCS observations can be attributed to aerosol effects. As circulation and vertical transport of dust are intertwined [Kahre et al., 2015], and presence of water ice clouds critically depends on temperature, complex feedbacks between modeled temperature and modeled fields of aerosols prevail, that need to be further explored in future model developments.    The improvement of mean temperature comes from the overall better representation of night and day temperatures, for most seasons. Two examples are shown in Figures 17  and 18, together with MCS data. In the equatorial region ( Figure 17) nighttime temperatures are reduced up to 10 K near the perihelion and about 5-6 K at NH autumn (L s =180 • ). However, daytime temperatures are warmer (∼ 8 K) than MCS around L s =250 • , and only a few K closer to MCS data for L s = 0-100 • . At high to middle northern latitudes (60 • N-80 • N) the LMD-MGCM reproduces MCS daytime zonal mean temperatures more accurately for all season, in particular in the NH autumn and winter, with data-model biases reduced up to 15 K from L s ∼150 • to L s ∼300 • (during dust season) while no significant differences are obtained in the equatorial region. Nighttime zonal mean temperatures at high to middle latitudes are closer to MCS observations for all seasons when the GW scheme is activated.

Improving the characterization of diurnal tides
Thermal tides are atmospheric waves forced by the diurnal cycle of incoming sunlight and they have a major impact on temperature and wind variability in the Martian atmosphere [Wilson and Hamilton, 1996]. The examples in Figures 19, 20, 21 and 22 show the quadrupole structure of the difference between dayside T pm and nightside T am temperature (T di f f = (T pm -T am )/2), as seen by MCS and predicted by the MGCM, for different solar longitudes, as indicated. The observed quadrupole of T di f f , centered roughly between 30 • S and 30 • N latitudes, is well known in the Martian atmosphere. It corresponds to the main response of the Hough mode of the diurnal tide, trapped between 22 • S and 22 • N latitudes and with a theoretical vertical wavelength of 30 km [Lee et al., 2009;Zurek, 1976]. This T di f f represents an important diagnostic quantity for the analysis of observations of the Martian atmosphere [Lee et al., 2009;Guzewich et al., 2012], which includes sun-synchronous tides and diurnal Kelvin waves. The primary effect of GW is to damp the thermal tides by reducing the diurnal oscillation of the meridional and zonal wind (see Section 5.2) Overall, the thermal tides were correctly represented by the model, but with differences in the values of the amplitude and vertical phasing when compared with MCS . The difference exceeded 15 K (the model being warmer than MCS between 1 and 10 −1 Pa), and the vertical phasing shifted slightly at higher altitudes. Figure 19, 20, 21 and 22 clearly show that the there are at least two seasons (e.g. NH spring and winter) during which the improvement associated with the inclusion of non-orographic GW parameterization is significant, especially the thermal tides amplitude above 10 Pa at the equator. The non-northern winter seasons see improvement in the depiction of the tide, but only right over the equator. It also seems that the wave depiction beyond the tide's latitude singularities is worse with the GW parameterization. The poorer performance in northern winter ( Figures  21 and 22) suggests that the changes to the zonal wind driven by GW may be spurious as the tide is ducted into the northern hemisphere by the strong zonal wind in this season, and the baseline GCM simulations seem to show this better, at least in the overall structure. With the non-orographic GW scheme activated, the maximum day-night difference value is around 10 K, comparable with observed values, and more interestingly the altitude of the peak of the tides is shifted down and its amplitude is also more realistic (between 22 • S and 22 • N latitudes, as observed). As discussed in the previous section, those improvements come from the overall better representation of both the nighttime (LT= 3:00) and daytime (LT= 15:00) thermal structure of the Martian atmosphere above 10 Pa.

Diurnal tide seasonal variation
The evolution of the main mode of the diurnal tide T di f f in the tropical band 20 • S-20 • N over all Martian seasons (MY29) is shown in Figure 23 for MCS data (top panel) and MGCM results (middle and bottom panels). The seasonal trend of T di f f is better represented by the MGCM with the GW scheme turned on (bottom panel), confirming the overall improvement in characterizing thermal tides in both amplitude and intensity, as described in the previous section. However, MCS data show a clear feature (e.g. the slight increase of tem-  Table 1 are used in our model simulations when the non-orographic GW scheme is on.  Table 1. perature around 20-30 Pa that shows up as a green blob) linked to a dust event at L s = 240 • which drastically changes the T di f f structure, and that the simulations presented in this paper fail to reproduce. This is probably related to the fact that the Martian thermal structure during dust storms is not well represented by the GCM. The storm is clearly visible in the Montabone et al. [2015] map of column dust optical depth for this MY, so perhaps the vertical dust distribution, size distribution, and/or water ice interactions are not correctly represented in our model, and they are producing the wrong thermal structure at that L s . Otherwise this suggests that during dust storms there are more GW propagating in the atmosphere, which could be taken into account in a future work by including a GW source varying with seasons in our scheme. It will be also crucial in the future to get the GW response to dust storms right when substantially larger storms are present (e.g. MYs 25, 28 and 34).

Sensitivity tests: impact on thermal tides
A number of sensitivity tests were performed to further evaluate the effect of GW induced drag on the GCM zonal wind and temperature and we use thermal tides (e.g. T di f f ) as a diagnostic. We found that our simulations, notably the zonal wind fields, were sensitive to all the set of parameters used in the scheme. Therefore the choice of the parameter ranges was done with caution, according to observational and theoretical constraints as described in Section 2.1. We select here 7 tests for a range of plausible values of gravity wave parameters and for different representations of the mean flow thermal forcing, listed in Table  2. For this set of runs the horizontal wavenumber k h range was fixed as in Table 1, corresponding to horizontal wavelength λ h between 10 km and 300 km (consistent with available observations); the other parameters vary with respect to the "best-fit" case in Table 1, which corresponds to Case 2 in Table 2). Specifically, the F 0 max was increased/decreased by one order of magnitude (Case 1 and Case 3 in Table 2, respectively) and the phase speed |c max | multiplied/divided by a factor 2 and 3 (Case 6 and Case 7, respectively). We found that the phase speed (c max ) is not a useful knob to tune compared with the momentum flux, which is instead one of the most important parameters in our scheme. For instance, values of F 0 max below 10 −8 kg m −1 s −2 produced a very negligible impact on the simulated field, while values larger than 10 −2 kg m −1 s −2 (not shown here) gave very unrealistic looking results. For details on the upper/lower limits of these parameters' probability distributions, see also Section 2.  Table 2 (see text for details). The contours are simulated day-night temperature differences T di f f between dayside (15 h) and nightside (3 h), shown for two seasons: Ls= 0 • -30 • N (upper panels) and Ls= 270 • -300 • (lower panels). The "best-case" simulation (Case 2 in Table 2), used as reference throughout this paper is represented in Panels C. MCS results for the same seasons (panels A) are also shown for comparison. Case 7 and Case 6, respectively. els B: Case 1) and increasing (Panels D: Case 3) of the probability range of the momentum flux amplitude at the source (F 0 max ) by one order of magnitude with respect the reference value (Panel C: Case 2). The larger the GW momentum deposition, the smaller the amplitude and the intensity of the thermal tide, and the more the main mode of diurnal tide is shifted at lower altitudes.
It is interesting to note that during the dust storm season Ls = 270 • -300 • (Figure 24 lower panels), stronger non-orographic GW momentum deposition (Panel D: Case 3) is required to better represent observations (e.g. an EP flux value between Case 2 and 3). This may indicate that winter jets add a further source of non-orographic GW and they could enhance the effect of waves propagation in the middle-upper atmosphere, suggesting that a source variable with season is required to fully reproduce the MCS fields (see also Section 6.2.1) A similar test was done varying the phase speed, more precisely varying the maximum value of the probability distribution of the propagating GW (c max ) from 10 m/s to 60 m/s, as shown in Figure 25 (Case 2, Case 6 and Case 7 in Table 2). Those values are comparable with background zonal wind ranges at the source level, and with typical values used in the literature [Medvedev et al., 2011]. As shown in Section 5, GW propagate vertically and break in the presence of strong zonal jets (the higher the contrast between the mean zonal flow wind and the phase speed, the higher the drag). In absence of wind jets, as in the southern hemisphere summer, higher phase speed (a factor 2 larger than the nominal values) reproduces slightly better the MCS results, since the waves can propagate vertically at higher altitudes before encountering critical levels and breaking. The best values identified for c max during these tests are between 30-40 m/s, which is within the realistic range. The last set of tests (Case 2, 4 and 5 as in Table 2) justifies the selection of the saturation parameter S c = 1 in our reference simulations (see Section 2). The saturation controls the altitude at which the GW break. Values between 1 and 5 (not shown here) produce similar results. Smaller values (e.g. S c = 0.1 as in Panel B) in our parameterization implies that only a small portion of wave energy is permitted to traverse the atmosphere above the source, and consequently it has a negligible impact on the GCM simulated fields (panels B can be compared with model results in Figures 19 and 21 for similar seasons.

Summary and conclusions
A stochastic parameterization of non-orographic GW was implemented into the Laboratoire de Météorologie Dynamique (LMD) Mars General Circulation Model (LMD-MGCM), following an innovative scheme, as described in Lott et al. [2012]; Lott and Guez [2013]. According to those authors, this approach is one of the best choices to represent the unpredictable aspects of the sub-grid dynamics, considering that convection and fronts can generate waves throughout the full range of phase speeds, wave frequency and vertical and horizontal scales. In addition, GW on Mars are still poorly constrained in terms of their basic parameters, including their sources of excitation, geographical distribution and vertical levels of dissipation. Given this uncertainty it is difficult to quantify the impact of non-orographic GW drag in the middle atmosphere of Mars using a unique set of parameters.
The LMD-MGCM already included a parametrization of orographic GW [Forget et al., 1999;Angelats i Coll et al., 2005] based on the GW drag scheme of Lott and Miller [1997] and Miller et al. [1989]. At the time the scheme parameters were chosen conservatively and the effect of the orographic wave drag probably underestimated. Here, we added a more general scheme which includes non-orographic GW, but it is possible that this also accounts for the previously underestimated effect of orographic GW. At this stage it is difficult to separate the contribution of orographic sources to the wind drag from non-orographic ones, but this will be certainly the next step in a future study. In our scheme the source is assumed to be located above typical convective cells (∼ 250 Pa) and a we have included a multiwave parameterization based on a stochastic approach.
The main goal here is to understand the role of non-orographic GW on the global circulation and the thermal structure of the Martian middle atmosphere. Our strategy was to identify a subset of "best-fit" GW characteristics by comparing LMD-MGCM simulations of the full MY29 and MCS/MRO temperature data, mainly between 1 and 0.01 Pa (∼ 50 and 80 km), where data-model biases were systematically larger before implementing the nonorographic GW scheme and locally reaching 10 to 20 K. The results show that the inclusion of a key physical mechanism such as the propagation of GW from the convective region to the upper atmosphere may partially explain those discrepancies at mesospheric layers, and it improves the accuracy of the LMD-MGCM, in comparison with MCS observations. Using the baseline parameters as in Table 1,the GW drag occurs mostly above 1 Pa (∼ 50 km altitude) and the main effect of GW is to slow down the winds everywhere. It predominantly reduces the westerly jets at high latitudes (negative drag) and weakens the equatorial easterly jet by accelerating the zonal wind (positive drag). Despite the GW decelerating the the westerly jets at mid-high latitudes (from tens to hundreds of m/s) in the NH winter solstice, with GW drag reaching maximum magnitudes of the order of 1 m/s/sol around 10 −2 Pa, the direct GW drag on the retrograde (easterly) wind at low latitudes is about two order of magnitudes weaker. This is due to the damping of the diurnal tides by the GW, and thus of the mean zonal wind, since the thermal tides -mean flow interaction is the main forcing of this retrograde jet.
Interestingly those results show that, in spite of the strong impact of GW on zonal winds, there are only minor changes on the predicted temperature field and a small reduction in the magnitude of the biases between MCS and LMD-MGCM at most seasons and latitudes. The physical explanation is that the changes in the temperature fields are not directly linked to the zonal mean drag but with the smaller GW-induced acceleration/deceleration of the meridional flow. There are however some notable exceptions of major changes in the temperature: (i) The mid-altitude negative bias at high southern latitudes in the first half of the year (northern spring and summer) is greatly reduced up to 20 K (ii) The mid-altitude positive bias at high northern latitudes in early winter (L s = 240-270) also decreases by 15-20 K (iii) Positive biases in northern mid-latitudes (40-60 • N) are generally reduced up to 10 K at all altitudes and times of the year, except at very high altitudes around local winter solstice (iv) Large negative biases at low altitudes at high latitudes (in northern fall) and northern mid-latitudes (in northern fall and winter) are not improved by the non-orographic GW parameterization and may be due to incorrect dust distributions and/or interactions between dust and water ice clouds during the dustiest time of the year. In addition, remaining biases with GW parameters may be attributed to aerosol effects, that need to be explored in further studies. Regarding the diurnal tides, they are better represented by the LMD-GCM when the GW scheme is on, especially their amplitude above 10 Pa at the equator. As previously discussed, the primary effect of GW is to damp the thermal tides by reducing the diurnal oscillation of the meridional and zonal wind. With the non-orographic GW scheme activated in the model the maximum day-night difference value is around 20 K, comparable with observed values, and also the peak of the tides is shifted at lower altitudes, in better agreement with MCS data. However, the simulations presented in this paper, (both with or without the non-orographic GW scheme included), fail to reproduce the slight increase of temperature linked to the dust event at L s ∼ 240 • observed by MCS around 20-30 Pa, which correspond to a global increase of day and night temperature. This may be related to the fact that the vertical dust distribution, size distribution and/or water ice interactions are not correctly represented in our model, and they are producing the wrong thermal structure at that L s . Otherwise this suggests that during dust storms there are more GW propagating in the atmosphere, and a non-orographic GW source varying with seasons may be required to further improve the results.