Investigating the Semiannual Oscillation on Mars using data assimilation

A Martian semiannual oscillation (SAO), similar to that in the Earth's tropical stratosphere, is evident in the Mars Analysis Correction Data Assimilation reanalysis dataset (MACDA) version 1.0, not only in the tropics, but also extending to higher latitudes. Unlike on Earth, the Martian SAO is found not always to reverse its zonal wind direction, but only manifests itself as a deceleration of the dominant wind at certain pressure levels and latitudes. Singular System Analysis (SSA) is further applied on the zonal-mean zonal wind in different latitude bands to reveal the characteristics of SAO phenomena at different latitudes. The second pair of principal components (PCs) is usually dominated by a SAO signal, though the SAO signal can be strong enough to manifest itself also in the first pair of PCs. An analysis of terms in the Transformed Eulerian Mean equation (TEM) is applied in the tropics to further elucidate the forcing processes driving the tendency of the zonal-mean zonal wind. The zonal-mean meridional advection is found to correlate strongly with the observed oscillations of zonal-mean zonal wind, and supplies the majority of the westward (retrograde) forcing in the SAO cycle. The forcing due to various non-zonal waves supplies forcing to the zonal-mean zonal wind that is nearly the opposite of the forcing due to meridional advection above ~3 Pa altitude, but it also partly supports the SAO between 40 Pa and 3 Pa. Some distinctive features occurring during the period of the Mars year (MY) 25 global-scale dust storm (GDS) are also notable in our diagnostic results with substantially stronger values of eastward and westward momentum in the second half of MY 25 and stronger forcing due to vertical advection, transient waves and thermal tides.


Introduction
In the upper stratosphere and mesosphere of Earth, the semiannual oscillation (SAO) of the mean zonal wind in the tropical stratosphere and mesosphere is commonly observed as a clear feature (Reed, 1966;Garcia et al., 1997). Reed (1966) used observational data to determine that the strongest westerly winds occurred shortly after the equinoxes in the lower mesosphere and then later extended to the lower levels. The pronounced SAO occurs at intermediate and upper levels of the stratosphere and mesosphere, and reaches a peak in amplitude at about the height of the stratopause.
Similar SAO-like features in the Martian tropics (between 10 • S -10 • N) have been noted in the work of Kuroda et al. (2008), mainly based on a free-running Martian Global Climate Model (MGCM).
Suspended dust particles are widespread in the Martian atmosphere, and they strongly affect its dynamical and thermal structure. Studying and characterizing the SAO phenomenon on Mars may contribute to understanding the dust cycle on Mars, as well as the Martian weather and climate.
So far, no comprehensive measurements of wind fields are available for the Martian atmosphere, although some localized snapshot measurements of middle atmosphere winds have been obtained from ground-based microwave measurements (Sonnabend et al., 2012). At middle and high latitudes, horizontal winds may be approximately computed using the thermal wind relation assuming geostrophic or gradient wind balance, via temperature fields observed by the instruments on board various satellites (e.g., Mars Global Surveyor (MGS) - Conrath et al., 2000, Smith, 2008and Mars Climate Sounder (MCS) -McCleese et al., 2010). However, even the gradient thermal wind relation becomes invalid at low latitudes so tropical winds are difficult to recover.
On the other hand, from variations in the structure and amplitude of solar thermal tides it may be possible to determine indirectly the corresponding variations of zonal-mean zonal wind. On Earth, Lindzen and Hong (1974) discovered that the zonal wind could affect the propagation of solar tides. Wilson and Hamilton (1996) also showed evidence of the effect of the solstitial mean zonal wind on the tide propagation on Mars. The seasonally varying wind should therefore also affect the pattern of solar tides. In the work of Kuroda et al. (2008), the SAO signal was evident by analyzing the solar thermal tides, derived directly from the observed temperature fields (from the Thermal Emission Spectrometer -TES onboard MGS) by taking the difference between the daytime and nighttime thermal fields. The thermal tides showed a clear SAO variation between 50 and 5 Pa in their study.
To conduct their detailed study of the SAO, Kuroda et al. (2008) made use of a simulation from a free-running MGCM. A clear SAO phenomenon between 10 • S -10 • N, similar in form to the stratospheric SAO on Earth, was described in their results (Figure 2 in Kuroda et al., 2008). To investigate the driving mechanism of the Martian SAO within the tropics, Kuroda et al. (2008) went on to calculate the different contributions to the tendency of the mean zonal wind, ∂u/∂t on interseasonal timescales. Their conclusion was that, at solstices, the horizontal cross-equatorial transport provided the westward forcing to the tropical zonal-mean zonal wind and vertical advection mostly provided a source of eastward momentum. The solar tides supplied eastward forcing at equinoxes. Somewhat surprisingly, at all seasons, Kelvin waves appeared to provide a westward forcing at all heights. Unlike on Earth, however, the effects of transient planetary waves and the eastward Kelvin waves on the Martian SAO were both quite weak.
Since the SAO may have further implications for the transport of tracers on Mars, as it has in the Earths stratosphere, it may be important to investigate the phenomenon of the Martian SAO in detail, with the additional aim of understanding the Martian dust cycle more completely. The intermittent nature of the available satellite measurements requires an interpolation in terms of time and space to acquire a complete observational coverage for the diagnosis of the Martian SAO. On the other hand, numerical models are able to overcome this difficulty to provide consistent data in both time and in three-dimensional space. A MGCM is therefore a perfect tool for analyzing the Martian SAO. The modeling study of Kuroda et al. (2008) seemed to provide a reasonably plausible representation of the Martian SAO in the tropics. But is a free-running MGCM, run with an unrealistically uniform dust distribution, sufficient to represent this aspect of the real Mars atmospheric environment? Various numerical modeling studies have proved capable of reproducing the full seasonal variability of the Martian climate, but MGCMs still struggle to fully represent the realistic interannual variability of the Martian climate system (Newman et al., 2002a,b;Basu et al., 2004Basu et al., , 2006Mulholland, 2012;Mulholland et al., 2013). In this context, data assimilation offers another approach towards obtaining a four-dimensional representation of atmospheric behavior that has the additional advantage of being consistent with both Martian observations and modeled physical constraints (Lewis and Read, 1995;Montabone et al., 2006;Lewis et al., 2007).
The goal of the present study is to diagnose the Martian SAO, based on the results from a Martian data assimilation system that is tightly constrained by real observations and can therefore provide a more realistic climatology than a free-running MGCM. The phenomenology and driving forces of the observed Martian SAO will be presented and investigated in this study. We wish to confirm if we observe a similar phenomenon in the tropics of the assimilated Martian climate measurements to that found in the work of Kuroda et al. (2008), and if the processes identified as driving the observed SAO are consistent with their analysis. A further question that can be answered by this observationallybased approach is whether components of an SAO can be observed outside of the Martian tropics.

Data and methodology
The dataset used in this study for diagnosing the Martian SAO is the MGS-TES Mars Analysis Correction Data Assimilation (MACDA) version 1.0 (Montabone et al., 2011(Montabone et al., , 2014. MACDA v1.0 is based on the UK-LMD MGCM that was first described in the work of Forget et al. (1999). A successive analysis correction scheme (Lorenc et al., 1991) was adapted and developed to assimilate retrieved temperature profiles and column-integrated dust opacities into the UK-LMD MGCM (Lewis et al., 2007).  Montabone et al. (2006) and Lewis et al. (2007) have already validated this dataset against available observations (other MGS-TES and radio occultation data), and this dataset has been proven to be reasonably consistent with those observations throughout this period.
MACDA is provided at a standard resolution equivalent to horizontal spherical harmonic spectral truncation T31, corresponding to a 3.75 • × 3.75 • nonaliasing dynamical grid. For this analysis, the  (Pike et al., 1984) is further used in each latitude band to reveal and decompose the phenomenology of the Martian SAO. In section 4, the Transformed Eulerian Mean equations (TEM) (Andrews and McIntyre, 1978) are applied to study the different contributions to the tendency of the zonally-averaged zonal wind in the tropical middle atmosphere throughout the seasonal cycle. The pressure and pseudo-height values used in the following figures were simply calculated from the model terrain-following vertical sigma coordinate, using a reference surface pressure of 610 Pa and a uniform assumed scale height of 10 km. It is noted that, because TES data only extends to ∼10 Pa (i.e. ∼40 km altitude), the reanalysis above that level is only indirectly influenced by observational increments being made below. Also, as the TES nadir data is only provided at the local time ∼02:00 and ∼14:00, this might not fully resolve the tides. Therefore, the reader should be cautioned that the tides in the reanalysis at altitudes above 10 Pa might mainly rely on the UK-LMD MGCM itself rather than on the data assimilated. In the Martian tropics (between 10 • N and 10 • S), the eastward-westward alternating pattern associated with the SAO was seen at altitudes between 100 Pa and 1 Pa, and most strongly at pressures above 30 Pa. This region is shown in more detail in Figure 2. The zonal wind begins with an eastward wind at the start of each MY, changing to a westward wind some time afterwards. This oscillation is evidently repeated twice a year, as in the experiments of Kuroda et al. (2008). In contrast to the cycles found in the uniform dust experiments by Kuroda et al. (2008), however, during each MY one of the two semi-annual cycles appeared to be stronger than the other. Clearly intensified eastward winds were also noticeable between L s ∼ 190 • and ∼230 • during the MY 25 global-scale dust storm (GDS). Normally, an eastward wind dominated below 20 Pa, but it was strongly reinforced and extended to higher altitudes (∼1 Pa) during this strong dust event.

Zonal wind
Apart from the tropics, the zonal wind within other latitude bands across the whole planet was also examined in the current study. It was clear from these that an SAO feature was also evident in latitude bands other than the tropics, though it was less pronounced towards the poles. The SAO pattern was also quite noticeable in the middle atmosphere (100 Pa -1 Pa) of the sub-tropics (10 • -40 • N and 10 • -40 • S) with the zonal-mean zonal wind direction oscillating from eastward to westward within a half-year cycle. Similarly to the tropics themselves, the two cycles in each year were of unequal magnitude, such that only one of the two cycles could be seen to extend throughout the whole atmosphere above 200 Pa. When the westward phase of the first SAO cycle dominated in the northern atmosphere, this corresponded to the westward phase of the second SAO cycle in the southern hemisphere extending vertically to the whole atmosphere above ∼200 Pa. This phenomenon was coincident with the annual cycle of the motion of the sub-solar point in latitude -which is located in the Northern Hemisphere in the first half-year, moving to the Southern Hemisphere in the second half-year. The sub-tropical SAO cycle in each hemisphere thus tends to be stronger when the sunlight is directly imposed upon that hemisphere. The SAO could also be seen in the middle and high latitude bands (60 • -90 • N, 40 • -60 • N, 40 • -60 • S, 60 • -90 • S) but much less strongly. This did not lead to a full reversal of the wind direction, as in the tropical atmosphere at altitudes higher than 30 Pa, but tended to slow down the prevailing eastward wind during the weaker cycle of the two in a given year.

Singular spectrum analysis
To further study the Martian SAO, a Singular Spectrum Analysis (SSA) was conducted on the daily-averaged zonal-mean zonal wind, in order to isolate the semiannual variations in each latitude band. SSA is a statistical technique applied in the time domain that is widely used in signal processing (Pike et al., 1984) and is similar to the Empirical Orthogonal Function (EOF) analysis which is usually applied in the spatial domain (e.g. Barnston and Livezey, 1987;Ghil and Mo, 1991) and decomposes the signal into separate contributions from a set of mutually orthogonal patterns in space and/or time.
Compared to other types of spectral analysis, the filters derived in SSA are not prescribed a priori, but are determined, optimally with respect to variance, from the data themselves (e.g. Collins et al., 1996) and therefore focus attention on signal components that contribute most strongly to the signal variance. Further, much more information can be obtained by studying the principal components of particular eigenvectors. SSA was employed here to demonstrate selectively the relative importance of any SAO signal compared to other oscillations, especially in those latitudes and altitudes where the SAO signal was superimposed upon other components.
In order to smooth the day-to-day fluctuations of the zonal wind fields, a 5-sol running-mean was applied to the diurnally-sampled zonal-mean zonal wind. Afterwards, the data was resampled in time by choosing every fifth data point for SSA analysis. The analysis time window here was chosen to SSA was carried out via the singular value decomposition of a covariance matrix, constructed by evaluating the time-averaged covariances of zonal mean zonal wind with respect to both altitude (pressure) levels and a set of discrete time lags, within each latitude band. The eigenvectors of such an SSA therefore represent two-dimensional recurrent patterns in both altitude (pressure) and time within the time window of the dataset. Eigenvectors are ordered with respect to decreasing eigenvalue, which is proportional to the corresponding with variance fraction "explained" by each eigenvector pattern.   Hemisphere exhibited signals with similar cycles to the SAO as in the third and fourth 310 PCs. It is interesting to note that the fifth and sixth PCs in the tropics (10°S -10°N) and 311

A TEM momentum budget for the SAO in the tropics
In order to identify the dominant processes generating the Martian SAO, the pressure-coordinate TEM equation for the zonal mean zonal wind (see e.g. Andrews and McIntyre, 1978) was applied to decompose the terms contributing to the acceleration of the zonal-mean flow. If all the terms contributing to the zonal-mean zonal wind tendency were moved to the right-hand side of the equation, this tendency can be written as follows, where u is the zonal-mean zonal wind, f is the Coriolis parameter, ϕ is latitude and r 0 is the radius of the planet, while the residual meridional circulation (0, v * , ω * ) is defined by The X in (1) represents unspecified zonal components of parameterized forcing or dissipation processes, such as friction or other nonconservative mechanical forcings.
The terms containing v * represent the acceleration of the mean zonal wind by meridional advection of momentum, while the second part (containing ω * ) represents the forcing due to vertical momentum advection. The third group of terms results from the non-zonal eddies, and is proportional to the divergence of the Eliassen-Palm flux (hereafter called the total wave forcing). These three contributions were computed respectively from the assimilated meteorological fields from MACDA and averaged in the tropics (between latitudes 10 • S and 10 • N).
To further investigate the roles of different zonal wavenumber components on the zonal-mean zonal wind, a two-dimensional (i.e. in the temporal and zonal dimensions, respectively) Fast Fourier Transform (FFT) was applied to the MACDA 2-hourly outputs to decompose the meteorological variables into various spatial and temporal harmonics. The contributions from different wave spectral components could be calculated, therefore, through the formulation of the EP flux terms using the decomposed meteorological variables. The different wave components decomposed in the present work consisted of 1) quasi-stationary waves with periods longer than 30 sols, 2) transient waves with periods longer than 1 sol and up to 30 sols, 3) westward thermal tides with periods of 1 sol or shorter, 4) eastward Kelvin waves with periods of 1 sol or shorter.
As the SAO phenomenon mainly happened above 100 Pa and up to altitudes of 1 Pa, the diagnosis of the forces driving it is focused on this part of the Martian atmosphere. Figure  In contrast, the vertical advection term (Figure 5b) mainly supplies eastward forcing to the zonal- The total forcing from different non-zonally symmetric waves (Figure 5c) exhibited more variation than the other two terms in (1). In MY 25 with a GDS, wave activity only provided westward forcing around 10 Pa in the first half of the year, while in MY 26 without a GDS, it also provided eastward forcing during L s ≈ 30 • ∼ 70 • . Apart from this, the total wave forcing exhibited a half-year cycle above ∼ 40 Pa. The phase of this oscillation was almost opposite to the forcing due to meridional advection above ∼ 3 Pa. Below 3 Pa altitude, the total wave forcing also partly contributed some westward forcing to the SAO of the zonal-mean zonal wind. The eastward forcing below 3 Pa started at L s ∼ 180 • , but it expanded throughout most of the second SAO. As a result, the eastward forcing due to total wave forcing appeared to support the eastward component of the SAO, but also weakened the westward forcing to some degree. Compared to the study by Kuroda et al. (2008), in the present analysis the eastward forcings in the second half of the year were obviously stronger in all three terms (meridional advection, vertical advection and total wave forcing). This is likely because of the increased diabatic forcing due to this being a more dusty season in the second half of the year (perihelion) than assumed by Kuroda et al. (2008) in their prescribed, smooth dust distribution. This effect was mentioned by Kuroda et al. (2008) as well, when they discussed the forcing due to meridional advection. Also, the eastward forcing in the vertical advection and total wave forcing terms reached their maximum around L s ∼ 220 • in the current study instead of ∼ 180 • in the study by Kuroda et al. (2008). This might be because the dustiest period in general does not usually happen at equinox, and this was not represented by a seasonally uniform dust distribution in a free-running simulation as was the case in the study by Kuroda et al. (2008). wind. Between 10 Pa and 2 Pa or below 50 Pa, the thermal tides provided eastward forcing. We note that they have also been found to lead to super-rotation within the height range of tidal forcing in the Martian atmosphere, given sufficiently dusty conditions (Fels and Lindzen, 1973;Lewis and Read, 2003). The eastward Kelvin waves (Figure 5g) consistently provided eastward momentum to u above 30 Pa. This is consistent with the study conducted by Moudden and Forbes (2008), in which the eastward propagating Kelvin waves were found to lead to strong eastward acceleration in the upper atmosphere associated with their dissipation. On the contrary, the westward forcing yielded by the eastward Kelvin waves in the study by Kuroda et al. (2008) is contradictory to the findings in this and other studies (e.g. Moudden and Forbes, 2008). A possible reason for this might be that the Kelvin waves in their model are forced non-physically, rather than generated spontaneously by the model physics. Below 30 Pa, the oscillation of the forcing due to eastward Kelvin waves in our analysis had a weak SAO signature, but with the opposite sign compared to the forcing due to meridional advection.
It is also worth mentioning the changes of the wave forcing during the GDS year, i.e. MY 25.
As generally the thermal tides dominated in the total wave forcing during the GDS, the significant increase of total wave forcing was mainly the result of increasing forcing due to amplified thermal tides.
Although the forcing due to quasi-stationary waves and transient waves exhibited some similarity in the annual pattern, they showed different responses to the GDS. The forcing due to transient waves supplied stronger eastward forcing during the GDS event, while the forcing due to quasi-stationary waves only led to a mild increase in the eastward forcing. However, in the year following the MY 25 GDS, the eastward forcing due to quasi-stationary waves was clearly stronger than it was in other MYs. In contrast, the GDS event had almost no impact on the eastward Kelvin waves above 100 Pa.
Finally, the residual of the actual zonal mean wind acceleration minus the sum of the zonal mean meridional advection, vertical advection and wave terms is shown in Figure 5h. This is smaller than each of the individual terms in Figures 5a-c below about 5 Pa (45 km), showing that the semi-annual pattern to the acceleration is largely explained by them. It is worth noting, however, that it is still significant in magnitude compared to the actual zonal wind accelerations (when averaged with a 60-sol window) and that in places, especially above 5 Pa, it is over an order of magnitude greater. This is likely to be a result of the strength of diabatic terms in the Martian atmosphere, in particular the short radiative timescale (about 1 sol) may lead to radiation acting as a damping in the residual term and making it more significant than in the case of the Earth. We also note the possibility that the larger residuals seen above 5 Pa may be, at least in part, an artifact of the free response of the model to the assimilation increments which cut off just below this level. Assimilation increments are only applied below about 40 km where TES retrievals exist, but we do not see any clear reduction in the size of the residuals when data is not assimilated and, in general, all the terms shown in Figure 5a-g are as large during these periods.

Disussion
In the current study, the phenomenology of the SAO in the zonal-mean zonal wind has been presented and reviewed. Unlike the preceding study of a Martian SAO by Kuroda et al. (2008), the present study used a dataset from a model directly constrained by observations to follow the actual weather during three Mars Years (Montabone et al., 2011(Montabone et al., , 2014Lewis et al., 2007). In the analysis, the SAO phenomenon appeared to happen not only in the tropics, but also extended its influence to higher latitudes. This oscillation appeared not always to manifest itself in completely reversing the direction of the zonal-mean zonal wind, but at some latitudes or pressure levels it was found to appear as the deceleration or acceleration of a prevailing unidirectional zonal-mean zonal wind. A Singular The diagnostic results here based on MACDA were somewhat different from the work of Kuroda et al. (2008), that was based on a free-running model simulation with a uniform assumed dust distribution. The observed dust distribution (Montabone et al., 2015) is clearly not consistent with such an assumption, so we should not be unduly surprised that this might lead to discrepancies between models and observations. The processes driving the tendency of zonal-mean zonal wind in the study of Kuroda et al. (2008) were mainly different from ours in the following ways: (1) the forcing due to the meridional advection in our study was less asymmetric between the two SAO cycles of a year, (2) the total wave forcing and decomposed thermal tides did not consistently display the oscillation as an SAO of zonal-mean zonal wind in every pressure level, (3) the strength of the forcing due to transient waves was comparable to that of the quasi-stationary waves, and (4) similar to the Earth, eastward Kelvin waves mainly supplied eastward momentum to the zonal-mean zonal wind, consistent with their dissipation of wave action. Finally, (5) the meridional advection and quasi-stationary wave forcing in the simulations of Kuroda et al. (2008) exhibited a peak or trough close to the surface which was not seen in the present work. The reason for this is not clear and should be investigated further in future work.
As the period of this dataset also included the period of the MY 25 GDS, the diagnosis also highlighted some distinctive features of the forces driving the circulation in this GDS year. The meridional advection supplied weaker westward momentum before the onset of MY 25 GDS, and the eastward momentum and westward momentum evident during the GDS were clearly stronger than those in other MYs. Along with the onset of the MY 25 GDS, the vertical advection and total wave forcing were also stronger. Among different wave components, the thermal tides were significantly intensified, most probably due to the enhanced solar heating effect due to increased dust during the GDS (see e.g. Lewis and Barker, 2005;Guzewich et al., 2014). Although the forcings due to quasistationary waves and transient waves in general displayed a similar pattern, the eastward forcing due to quasi-stationary waves appeared to be stronger in the year after the MY 25 GDS instead of during the MY 25 GDS itself.
In future work, the TEM analysis should be extended to the latitude bands outside of the tropics, making it possible to discover differences among the different processes across different latitudes. To establish any further link between the features of different forcings and the trigger/response of a GDS, it would be desirable to investigate at least one other GDS event. This will require the extension of reanalysis datasets, such as MACDA, to include MYs with other GDSs, such as occurred in MY 28.

Data availability
The MACDA dataset on which this paper is based can be accessed via the website of the Centre for Environmental Data Analysis of the UK Natural Environment Research Council (https://catalogue.ceda. ac.uk/uuid/01c44fb05fbd6e428efbd57969a11177) and is documented by Montabone et al. (2014).