Atmospheric thermal tides are planetary-scale gravity waves, with periods that are harmonics of the solar day (Chapman and Lindzen, 1970; Forbes, 1995). Tides are particularly prominent in the Martian atmosphere so that temperature, wind fields and surface pressure exhibit strong dependence on local solar time (LT). This is because the atmospheric density near the Mars surface is very low, two orders of magnitude smaller than on Earth, and the Mars surface thermal inertia is low compared to that of the terrestrial oceans. These factors lead to very large diurnal variations in surface temperature and a resulting strong atmospheric response. The rapid radiative damping time results in a strong coupling between surface and boundary layer temperature. In addition, aerosols also lead to significant solar heating within the atmosphere, particularly during dust storm episodes (Zurek, 1976; Zurek et al., 1992; Wilson and Hamilton, 1996).
Thermal tides are particularly significant for the understanding the Martian atmosphere due to their broad impacts on Martian atmosphere dynamics (Forbes et al., 2002), and they also have impacts on dust and water transport, which can in turn affect the atmospheric stability, dynamics, chemistry and thermal structure (Zurek et al., 1992). Tides are also of great importance for interpreting spacecraft observations, especially for spacecraft in sun-synchronous orbits.
Tides can be decomposed into sun-synchronous westward propagating, i.e. migrating waves that are strongly driven by solar heating, and stationary, i.e. non-migrating waves resulting from zonal variations in thermo-tidal forcing. Zonal modulation of forcing can arise from longitudinal variations of topography (the most significant factor for zonal modulation of tides), albedo and surface thermal inertia, as well as from spatial variations in radiatively active aerosols including dust and water ice clouds (Zurek, 1976; Wilson and Hamilton, 1996).
The study of Martian thermal tides and stationary waves has commonly been based on temperature observations from the Mars Global Surveyor (MGS) Thermal Emission Spectrometer (TES), which provides full latitude and longitude coverage, as well as in situ measurements at the Viking (and now Curiosity) landers. TES observes mostly in nadir viewing mode in a sun-synchronous orbit that provides twice-daily temperature observations (2 am and 2 pm local solar time, LT), resulting in atmospheric soundings from 0 to 35 km (surface to 10 Pa) with a vertical resolution of ~10 km (Smith et al., 2001). More recently, the Mars Reconnaissance Orbiter (MRO) Mars Climate Sounder (MCS) (McCleese et al., 2007) has provided considerable additional insights to Martian thermal structure (e.g. Lee et al., 2009; Guzewich et al., 2012; Kleinböhl et al., 2013); in this study we focus on TES, reserving results using MCS for a future publication.
The twice-daily TES observations at 2 am and 2 pm LT provide some constraint on the tides, but are insufficient to fully characterise the basic features of the thermal tides. Instead, for Mars tide studies, a data assimilation system optimally combines temporally and spatially irregular observations with a background state from a global climate model (GCM) forecast. Mars data assimilation systems have already been developed to yield a dynamically consistent set of dynamical and thermal fields providing the basis for detailed investigations of the Martian dynamical features and circulations system (Houben, 1999; Zhang et al., 2001; Montabone et al., 2005; Lewis et al., 2007; Hoffman et al., 2010; Lee et al., 2011; Greybush et al., 2012; Navarro et al., 2014; Steele et al., 2014). The resulting analyses give an estimation of the evolving state of the Martian atmosphere during the MGS mission and provide a basis for analysing Martian atmospheric variability, circulation features, dust transport, etc.
Lewis et al. (2007) performed the first complete TES-based reanalysis, which made many contributions to the understanding of Martian atmospheric properties such as atmospheric tides (Lewis and Barker, 2005), dust storm variability (Montabone et al., 2005), tropical water ice clouds (Wilson et al., 2008) and predictability (Rogberg et al., 2010). Since then, ensemble-based methods have developed and matured (see Kalnay and Yang, 2010), and have been identified as a key advancement in meteorological data assimilation. Ensemble-based methods estimate the flow-dependent background error covariances, i.e. the multidimensional background uncertainty, from the ensemble spread, i.e. variability, and this has led to significant improvements in Numerical Weather Prediction (NWP) skill (e.g. Kleist, 2012). Ensemble data analysis exploits inter-variable correlations to estimate both observed and unobserved model variables (e.g. updating surface pressure and winds from observations of temperature). An example of an ensemble-based data assimilation system used in Martian atmosphere reanalysis is the Data Assimilation Research Testbed (DART) (Anderson et al., 2009) used by Lee et al. (2011), where TES radiances were assimilated into the MarsWRF GCM (Richardson et al., 2007). The depiction of some atmospheric features such as the circulation, the CO2 cycle and surface properties were shown to be improved through this ensemble-based assimilation. Montabone et al. (2014) assimilated MGS TES temperature retrievals and column-integrated dust opacities into the UK-LMD-MGCM in use at the University of Oxford and at the Open University in the UK. Following Lewis et al. (2007), this assimilation scheme uses the Analysis Correction (AC) scheme (Lorenc et al., 1991), and the successive corrections method that demonstrated its robustness in Mars reanalysis, even during the imperfect MGS aerobraking period. A key difference between the AC scheme and others (such as variational and ensemble methods) is that the former nudges a model forecast gradually towards observations rather than making a correction towards observations at discrete time steps. The resulting Mars Analysis Correction Data Assimilation (MACDA) provides a multi-year reanalysis spanning almost three complete Martian seasonal cycles (Montabone et al., 2012). An advanced version assimilating both MGS TES and MRO MCS (McCleese et al., 2007) profiles into the UK-LMD-MGCM provided the first results of a complete assimilation of both datasets (Montabone et al., 2014), covering almost six complete Mars years (MY), and allowing for the study of the interannual variability of the Martian weather and climate with respect to its components including the thermal tides.
Hunt et al. (2007) developed the local ensemble transform Kalman filter (LETKF) system, an advanced data assimilation system that has been successfully tested in the terrestrial data assimilation (Miyoshi and Yamane, 2007; Szunyogh et al., 2008; Miyoshi et al., 2010). The LETKF optimally combines the ensemble forecasts in order to optimally fit the observations. It has several advanced features such as estimating observation error covariance inflation parameters (Li et al., 2009), a no-cost smoother (Yang et al., 2012) and an advanced quality control based on Ensemble Forecast Sensitivity to Observations (Hotta, 2014). Greybush et al. (2012) used a 4D-LETKF scheme that assimilates observations at their time of observation, rather than in 6-hr windows centred at the assimilation time (as had been used in Earth studies), coupled with the Geophysical Fluid Dynamics Laboratory (GFDL) Mars Global Climate Model (MGCM; Wilson, 2011), and generated a new analysis of Martian weather and climate. The 6-hr forecast root mean square difference (RMSD) with respect to observations were found to be much smaller in the 4D-LETKF assimilation than in a free running MGCM, suggesting a successful assimilation. However, we subsequently found that the assimilation resulted in eastward-propagating semi-diurnal tides with unrealistically large amplitudes. One of the main results of this paper is that these unrealistically large tides are due to a spurious resonant forcing associated with the use of 6-hr assimilation windows. We then demonstrate a solution that allows the creation of a Mars reanalysis with realistic diurnal features.
In this study, the LETKF coupled with the MGCM is used to assimilate TES observations. Different versions of the data assimilation system with different data windows, different parameterisations of clouds etc., yield different estimates of the features of the thermal tides. Background information on the TES observations, Viking Lander and Curiosity data, the MGCM and LETKF, and a basic description of tides are presented in Section 2, the configuration of the system with different assimilation windows is discussed in Section 3, results are described in Section 4, and Section 5 is a summary and discussion.
The MGS spacecraft was launched in November 1996, reached Martian orbit in September 1997, started its primary mission in February 1999 and failed in November 2006 during its third extended mission. The MGS spacecraft followed a sun-synchronous orbit during its primary science mission, providing 12 daily observations each at 2 am and 2 pm local solar time (LT), drifting ~30° in longitude with each subsequent orbit. In the current study, we assimilate temperature retrievals into the MGCM, obtained from TES radiance measurements from 1997 to 2004 (MY 24 to early MY 27) (Christensen et al., 1992; Christensen et al., 2001). Atmospheric temperature soundings extending from the surface to ~40 km (0.1 hPa) with a vertical resolution ranging from ~10 km near the surface to ~15 km at higher altitudes are retrieved from nadir radiances in the 15 micron CO2 absorption band (Conrath et al., 2000). The TES nadir radiances are also used to retrieve column dust and ice opacity, but these are not used in this study.
TES temperature errors are estimated to be around 3 K, with higher errors in polar regions and near the surface (Eluszkiewicz et al., 2008). Systematic errors can also vary from orbit-to-orbit. During dust storm seasons, retrieval uncertainties of both temperature profile and dust properties are relatively high, partially due to dust being assumed to be composed of non-scattering particles (Smith et al., 2001).
The Viking program sent two space probes to Mars, each consisting of an orbiter that obtained high-resolution images of the Martian surface and a lander that recorded the Martian atmosphere data from the surface. Viking 1 was launched on August 20, 1975, entered Mars orbit on June 19, 1976. The Viking Lander 1 (VL1), landed on July 20, 1976 at 22.48°N 49.97°W and retired on November 13, 1982. Viking 2 was subsequently launched on September 9, 1975, entered Mars orbit on August 7, 1976. Viking Lander 2 (VL2) landed on September 3, 1976 on 47.97°N 225.74°W and retired on April 11, 1980 (Kieffer et al., 1992). These two Viking Lander sites provided surface pressure time series consisting of an almost 4 Martian year cycle, allowing the study of the seasonal variation of Martian climate. Records of thermal tides at these two Viking Lander sides show a lack of interannual variability in the aphelion (furthest from the Sun) season, which is consistent with the high degree of repeatability in this season for atmospheric temperature observed by Viking and MGS (Liu et al., 2003; Smith, 2004). In this season, radiative processes force the Mars atmosphere and the dynamics are non-chaotic. Thus the repeatability of Martian tides in this particular season may be attributed to a repeatable variation of atmospheric aerosol and therefore repeatable thermo-tidal forcing as modulated by the topography and surface thermal inertia. The small interannual variability supports the comparison between the Viking data and the TES assimilation in the present study. A harmonic analysis of the surface pressure tidal components for the two Viking Lander sites has been conducted by Wilson and Hamilton (1996) and Bridger and Murphy (1998), suggesting both phase and amplitude of the diurnal and semi-diurnal surface pressure tides at the two Viking Lander sites have systematic seasonal variations, which are consistent with the presence of resonantly enhanced Kelvin waves due to the dynamical effects induced by topography forcing (Wilson and Hamilton, 1996; Bridger and Murphy, 1998). A detailed discussion regarding this issue is provided in Section 4.1 and 4.3.
After detaching from the Mars Science Laboratory (MSL), the Curiosity Rover landed within the Gale Crater (4.59°S 137.44°E) on August 6, 2012. Curiosity's Rover Environmental Monitoring Station (REMS, described in detail at www.msl-scicorner.jpl.nasa.gov/Instruments/REMS/ and by Gómez-Elvira et al., 2012) has measured pressure, temperature, wind, humidity and ultraviolet fluxes since that time. Hourly data from REMS is used here as an independent verification. Note that the rover sits within the Gale Crater at a distinctly higher surface pressure than is represented by the model results, which will have some influence on the comparison of simulations with observations (Haberle et al., 2014). As discussed in Tyler and Barnes (2013), the crater appears to enhance the tide response, and this matter is discussed in Section 4.3.
The MGCM was originally based on the GFDL SKYHI terrestrial GCM (Wilson and Hamilton, 1996) and is now adapted to the GFDL Flexible Modeling System (FMS), using a finite volume dynamical core (Lin, 2004). The MGCM has already been used in Mars atmospheric tides and planetary wave analysis (Wilson and Hamilton, 1996; Hinson and Wilson, 2002; Wilson, 2002; Kleinböhl et al., 2013). The MGCM includes surface and subsurface physics essential for the surface temperature calculation; a budget containing gaseous and condensed CO2 that tracks the evolution of global average surface pressure through the annual cycle; inventories of water (regolith, surface frost, water vapour and condensate) and dust aerosol that maintain mass conservation, and a Mellor-Yamada 2.5 boundary layer scheme. Wilson (2011) introduced a parameterisation of subgrid-scale gravity waves forced by topography.
In this study, we use essentially the same baseline model configurations as used previously with the LETKF by Hoffman et al. (2010) and Greybush et al. (2012). This configuration is 5×6° (or ~300 km) resolution with a 36×60 latitude-longitude horizontal grid. A hybrid pressure-sigma vertical coordinate system is employed in this model, with 28 vertical levels. A higher resolution (2×2.4°, or ~120 km) simulation was also used in Greybush et al. (2012) to test the model sensitivity to horizontal resolution: the features found at both resolutions in most of the free atmosphere such as the polar westerly and tropical easterly jets, and the zonal mean temperature structure are very similar. Although the model is mass-conserving, data assimilation can lead to spurious (i.e. unphysical) sources and sinks. To maintain CO2 mass conservation in the assimilations, at each grid point the analysis surface pressure is scaled by the ratio of background atmospheric mass to analysis atmospheric mass after each analysis step, which allows local surface pressure to match the temperature increments while conserving atmospheric CO2 (Greybush et al., 2012).
Dust aerosol is a major source of atmospheric heating of the Mars atmosphere. A ‘seasonal dust’ configuration in which the dust mixing ratio vertical profile is specified as in Forget et al. (2001) and Montmessin et al. (2004), is used in several assimilation experiments in this study. Simulations with an interactively evolving three-dimensional dust distribution were also carried out, with the opacity field represented by radiatively active dust tracers that include advection and sedimentation. However, such an evolving three-dimensional dust distribution did not yield qualitative differences in the tidal analysis from the analysis using ‘seasonal dust’ during the time of year of our study, therefore to simplify our experiment we only use ‘seasonal dust’.
Water ice clouds also make a significant radiative contribution in the Mars atmosphere (Wilson et al., 2008; Madeleine et al., 2012). The water ice cloud microphysics scheme employed in this model is as described in Montmessin et al. (2004). Cloud microphysical processes are critical in this model. The microphysical processes consider transportation of cloud ice mass and predict a spatial mean cloud particle size by specifying the ice mass and the number (N) and size (ro) of dust-like cloud condensation nuclei. Assumptions associated with the background dust particle size distribution (PSD) impact the cloud particle size, which controls the sedimentation and cloud optical properties and eventually determines cloud vertical distribution. Another parameter involved in this study is a scaling Constant for Radiation (CR) that determines the cloud opacity. The introduction of CR facilitates the exploration of cloud radiative impact due to the straightforward switch between simulations with radiatively passive (CR=0) and radiatively active clouds (CR>0). As discussed in Greybush et al. (2012) and Kleinböhl et al. (2013), we can vary CR to ‘tune’ the cloud radiative impact.
Data assimilation methods optimally combine information from current observations and a prior short-term forecast, creating an analysis that represents the best estimate of the state of the atmosphere. This state estimate is then used to initialise the next short-term forecast, which then serves as the ‘first guess’ or the background, for the next analysis. In the current study, the forecast model is the MGCM. In the analysis, the background error covariance (denoted B) and the observation error covariance (denoted R), characterising the background and the observation uncertainties determine the optimal weighting of model forecast and observations. The background error covariance (B) not only quantifies the forecast uncertainty but also determines how one type of observation influences other model forecast variables through the covariance between different variables. In this study, only temperature observations are utilised to update other model variables such as surface pressure and wind. Hoffman et al. (2010) showed the effectiveness of assimilating temperature retrievals in reducing the RMS difference with observations compared to the freely running model in an Observing System Simulation Experiments (OSSE). As in Hoffman et al. (2010) we adopted the LETKF scheme (Hunt et al., 2007). In our experiments the ensemble members differ both in terms of initial conditions and dust and water ice cloud properties.
While obtaining comparable performance as other sequential ensemble-based Kalman filters applied using the same NWP model (Whitaker et al., 2008), the LETKF is computationally efficient because the LETKF assimilation is performed independently at each grid point in the centre of a local region, making the LETKF computations completely parallel.
The LETKF algorithm is now briefly described (see Hunt et al., 2007, for a detailed derivation) In the background ensemble of model states, the kth member of the total K ensemble members is represented as a vector and the ensemble mean is denoted as . Then the column k of Xb denotes the difference between model state of the kth ensemble member and the ensemble mean :
Then the background error covariance matrix B is expressed as
Similarly, the local analysis error covariance is written as
where Xa represents the analysis ensemble perturbation matrix, and , the analysis error covariance in ensemble space (dimension K×K), is given by
where Yb refers to the background ensemble perturbation matrix transformed into the observation space with the observation operator H; R is the observation error covariance, assumed to be diagonal for simplicity (i.e. the observation error covariance between different variables is neglected). I is the identity matrix of rank K. Generally, the ensemble size K is much smaller than the model dimension, and as a result, the calculation in ensemble space greatly enhances the computational efficiency.
Since we have a limited ensemble size, sampling errors due to spurious long-distance correlations are inevitable. Localisation is used to reduce this kind of sampling error, as well as to reduce the computational cost. Within the LETKF scheme, ‘R-localisation’ is achieved by multiplying the observation error covariance matrix R by fR and ignoring observations beyond a cutoff distance, which, in this study, is chosen to be 3.65 times the horizontal localisation distance L following Miyoshi and Yamane (2007). Here fR is specified by an inverse Gaussian function that exponentially decreases the observation impact up to the cutoff distance:
Here L is the localisation length scale, and dij the distance between the ith observation and the jth model grid point (Hunt et al., 2007). In this way, the uncertainty assigned to the observations (observation error) increases smoothly until they have very large uncertainty and therefore no influence on the analysis. Greybush et al. (2011) demonstrated in examining the wind-height geostrophic balance that the performance in terms of analysis accuracy of the LETKF using R-localisation is similar to other ensemble Kalman filters using B-localisation (where localisation is applied to the background covariance matrix B) (Hamill et al., 2001; Houtekamer and Mitchell, 2001). The choice of localisation scale L is one of the topics that we investigate, and is discussed in Section 4.1. The analysis ensemble mean is obtained as
where the weight for the ensemble mean is given by
The analysis ensemble perturbations are obtained from:
Let be the columns of the matrix that are the result of adding the ensemble mean to each of the columns of Wa, The resulting ‘weight’ vectors then give the analysis perturbations as linear combinations of the background ensemble perturbations, so that:
The linear combination of the background ensemble perturbations can be seen as the analysis increments, in which the relative contribution of each ensemble member to the analysis is available. If each ensemble member contains a unique configuration of model parameters [in this study, the optimal dust scale and the radiatively active water ice cloud (RAC) content], eq. (9) can be used to assess the most appropriate physical parameter that best fits the observations. Since the model deficiencies and sampling errors are frequently underestimated in the ensemble Kalman filters, the background ensemble spread should be artificially increased to allow the LETKF to give more weight to the observations due to the probable underestimation of the uncertainties of the ensembles. This process, called ‘inflation’, is discussed in Section 3.
As indicated in the introduction, tides are a dominant aspect of Mars meteorology and are much more prominent than those on Earth. Unlike Earth tides, Mars tides are strongly influenced by aerosols (dust and water ice clouds). The longitude-time dependence of thermal tides can be represented as:
and can be expressed in a local time reference frame as:
where s is the zonal wavenumber, λ denotes east longitude, σ represents the temporal harmonic (σ=0 refers to a stationary wave, not discussed in this study, σ=1 stands for the diurnal tide; σ=2 represents the semi-diurnal tide), tLT means the local solar time and δs,σ is the phase. Universal and local solar time are related by t=tLT−λ. Tides with s>0 (s<0) correspond to westward (eastward) propagation, and s=0 means zonally symmetric tides. When σ=s, tides are sun-synchronous without longitude dependence in the local time coordinate system and are also called ‘migrating’ tides. Tides that do not follow the motion of the sun (σ≠s) are called ‘non-migrating’ tides, which typically result from thermo-tidal forcing zonal variations. In a local solar time reference, an observed wave number m zonal wave may represent a mix of stationary Am,0 and non-migrating tide components [As,σ, with (s−σ)=±m] (Forbes and Hagan, 2000; Wilson, 2000). For example, an observed wave 2 variation includes contributions from the wavenumber 3, westward-propagating diurnal tide (A3,1) and wavenumber 1, eastward diurnal tide (A-1,1), as well as higher temporal harmonics (A0,2, A4,2, A1,3, A5,3, …).
Classical tide theory can be used to understand the latitude and vertical structure of the atmopsheric response to imposed thermal forcing. The diurnal response can be represented as a series of vertically propagating waves in the tropical region and vertically trapped waves in the extratropics (Chapman and Lindzen, 1970). The dominant tropical diurnal tide response is a vertically propagating component with a ~33 km vertical wavelength (Zurek, 1976; Wilson and Richardson, 2000), while other westward non-migrating tides consist of a series of shorter vertical wavelengths. Since the vertical wavelength of diurnal tide in tropical region and the vertical smoothing length of the TES nadir retrievals have the same order of magnitude, the observed tidal amplitudes are significantly reduced compared to the atmospheric response (Wilson and Richardson, 2000). By contrast, the diurnal tide response in the extratropical region is vertically trapped and so varies little over the depth of the forcing region.
There is a class of Kelvin waves that is particularly prominent on Mars. As noted in Hamilton and Garcia (1986) and Zurek (1988), the external mode on Mars is such that the free mode period for Kelvin waves is roughly 24 hr, rather than 33 as on Earth. The proximity of this period to the solar day is what makes resonant enhancement much more likely on Mars than on Earth. The Kelvin mode contributions to A-s,s (s=1,2,3 …) are relatively non-dispersive and are thereore all potentially nearly resonance for forcing at 24/s periods.1 The detection of the near-resonant Kelvin waves in the atmosphere including a comparison with a MGCM simulation has already been discussed in previous study (Hinson et al., 2008). They can be readily excited by the scattering of the migrating tide by planetary-scale topography. Attention has been paid to the most prominent components of the eastward-propagating waves, in particular to diurnal Kelvin waves (DK1, DK2, … corresponding to s=−1, −2, …) which are solutions of the Laplace Tidal Equation. They show meridionally broad symmetric structures (Chapman and Lindzen, 1970). Strictly speaking, A-s,s can be represented by a number of different Hough modes with different meridional structure. The Kelvin wave is the dominant mode with the broadest meridional structure. The vertical structure of DK1 resembles well the equivalent barotropic Lamb wave and may be enhanced by resonance (Zurek, 1976, 1988; Wilson and Hamilton, 1996). Both DK2 and DK3 are vertically propagating waves with ~80 and ~60 km vertical wavelengths, respectively; the amplitudes of both grow exponentially with height (Wilson, 2000). The zonal wavenumber 2, semi-diurnal east wave (SE2), labelled as A-2,2, could be resonantly enhanced due to different factors, i.e. assimilation window length, and is a major point discussed in this study. We are also interested in the westward-propagating, wavenumber 2 semi-diurnal tides (A2,2), the dominant contribution to semi-diurnal wave 2 (S2). This tide has a very long vertical wavelength and a broad meridional structure, so it effectively integrates tide forcing distributed over a large range of heights and latitudes (Wilson and Hamilton, 1996; Forbes and Miyahara, 2006).
A difficulty in assessing the accuracy of the tides in the analysis fields is the sparseness of observations. A fundamental issue for data assimilation is whether the twice-daily spacecraft observations are sufficient to constrain the tides in the assimilating model. One hypothesis is that that the migrating tides are not really constrained by the observations though the non-migrating tides may be more readily constrained. In the current study, the time series of temperature and surface pressure are decomposed into a sequence of harmonics of a solar day. As in the study of Forbes and Miyahara (2006), in which the fully resolved diurnally varying surface pressure observations provided at the Viking Lander and Curiosity sites are used to study the semi-diurnal tides, the tidal signals in the analyses are evaluated at the Viking Lander and Curiosity sites, and are also compared to the observations to get quantitative assessments (Section 4.3).
Experiments assimilating TES profiles are conducted for a 20-sol period during the Northern Hemisphere (NH) Martian summer of MY25, approximately Ls (solar longitude) 94.5°–103.7°, or sols 388–408 past the perihelion. The TES data were acquired from the Planetary Data System (PDS). Each year during the NH Martian summer (aphelion), atmospheric variability is primarily governed by the thermal tides and aerosol structure, as travelling wave activity is at a minimum and dust opacities are relatively small (0.15–0.3) in most of the atmosphere. Allowing for a 9-sol spin-up period, results are compiled for sols 397–408 (Ls 98.6°–103.7°).
Greybush et al. (2012) conducted experiments with both 16 and 32 ensemble members, but found the improvement in forecast accuracy with 32 members not significant. As a result, we use 16 ensemble members in this study to save computational costs. A dust configuration with a prescribed vertical profile of dust mixing ratio that evolves with the season and latitude according to an analytic formula (Montmessin et al., 2004), which is labelled as ‘seasonal dust’, described in Section 2.2, is applied in the current study. Using a dust scenario that is uniform in longitude makes it easier to diagnose temperature changes resulting from longitudinally-varying analysis increments from the data assimilation. In this study, the dust amplitude increases linearly and the dustiest ensemble member has 2.5 times more dust than the least dusty ensemble member.
Observations that are geographically close enough are combined into ‘super-observations’ (Alpert and Kumar, 2007). This significantly reduces the number of observed temperature profiles assimilated, reduces computational time, filters out small scale features that can not be represented by the model, and reduces random observation error (Greybush et al., 2012). Only observations with good quality control flags are assimilated. Quality control is also applied within the LETKF, rejecting observations if the departure from the background is larger than a QC threshold. We use the same configuration in the current study as in Greybush et al. (2012): the prescribed TES observation error is 3 K, and the QC threshold is set as five times the observational error.
Like other EnKF systems, the LETKF system is overconfident about the quality of the ensemble forecast because it neglects parts of the model errors and the impact of nonlinearities and thus underestimates the analysis and forecast ensemble spread. A significant source of model error is the poor knowledge of radiatively active aerosol in the atmosphere, including dust and water ice clouds. This reduction of ensemble spread results in an underestimation of the weights of the observations and requires ‘inflating’ the spread at each assimilation time. We use adaptive inflation (Miyoshi, 2011) in our assimilation, in which a Kalman filter temporally smooths the innovation statistics (i.e. O-B statistics) to estimate multiplicative inflation values at each grid point (Hunt et al., 2007). Since the satellite path only covers a small fraction of the whole domain, we only inflated the ensemble analysis at the grid points where there are observations – at other locations with no nearby observations, the ensemble spread is not reduced during the analysis step. To maintain dynamic consistency, the inflation is applied to perturbations of zonal and meridional wind components, as well as temperature, but not to surface pressure, because that would not conserve the mass of the atmosphere. For TES profiles, as no observation information available above ~0.l hpa (10 Pa) to constrain the ensemble spread, inflation is not needed in these levels, and a linear decrease in spread inflation values from the top level of the TES observations is applied. Since we use different assimilation window lengths, the initial adaptive inflation values are scaled proportional to the assimilation window length. In this study, for 6-hr LETKF assimilating TES profiles, covariance inflation value is initially 77% at model level 7 (the top level of TES observation) and decreases linearly to 0% at model level 4 and the levels above where there are no observations. Since the background spread is inflated three times within 6 hr for 2-hr window assimilation, covariance inflation value is adjusted to 21% . Similarly, covariance inflation value is adjusted to 10% for 1-hr assimilation window.
In this study, we begin with the same model configuration as in Greybush et al. (2012): horizontally 600 km localisation, vertically 0.4 log P localisation and temporally 3-hr localisation. We have found that the semi-diurnal east wave response can be sensitive to the horizontal localisation radius, which will be discussed in Section 4.1. Post-processing steps applied at each analysis step include a Shapiro filter in space to remove high frequency noise in the analysis increments at high latitudes. As with terrestrial data assimilation systems, a filtering to the analysis increments is needed to help reduce the generation of spurious gravity waves in forecasts launched from the analyses. Greybush et al. (2011) discusses the impact of localisation on balance in EnKFs, and that imbalance can be minimised (but not eliminated) by selecting a localisation length scale that is sufficiently long. As Mars and Earth have similar Rossby radii of deformation, similar reasoning should apply on the red planet, and our 600 km half-length localisation length scale is similar to that used for the SPEEDY model on Earth (with a comparable horizontal resolution). In addition, enforcement of CO2 mass conservation is obtained by scaling the analysis surface pressure at each grid point by the ratio of the background atmospheric mass to the analysis atmospheric mass.
Figure 1 shows the spatial patterns of the semi-diurnal surface pressure amplitude (in percent) for a free run simulation and data assimilation with a 6-hr window. In all cases, the surface pressure field is first normalised by the diurnal mean surface pressure to account for topography (in the spirit of converting to ‘sea level pressure’ on Earth). An examination of the tide components confirms that the semi-diurnal tide response is dominated by two modes, the migrating tide (following the Sun) and semi-diurnal east wave (SE2) (Wilson and Hamilton, 1996), with their amplitudes as a function of latitude shown in the bottom row of Fig. 1. Wave number 4 patterns (Fig. 1a and b) can be seen as the interference between these two dominant semi-diurnal modes (Fig. 1c and d), the westward- and eastward-propagating tides, A2,2 (the migrating tide) and A-2,2 (the SE2 tide) as described in Section 2.5. Compared to the free run spatial patterns, the 6-hr TES assimilation shows stronger wavenumber 4 patterns. At Ls~100, the surface pressure oscillation at Viking Lander 1 (VL1, 22.7°N 48.2°W) has an amplitude of about 0.3–0.4% (Wilson and Hamilton, 1996). The 6-hr TES assimilation shows a much too strong signal (2.5%) compared to the Viking Lander observation data. Meanwhile, the surface pressure semi-diurnal tide (S2) amplitude at Curiosity site (4.59°S 137.44°E) has a calculated value of 0.5–0.6% by tidal fitting, and the 6-hr TES assimilation again shows a much stronger signal (1.9%). The enhanced wave number 4 patterns in the TES assimilation (Fig. 1b) is due to the significantly enhanced amplitude of A-2,2 (SE2) (Fig. 1d).
It is clear that the 6-hr assimilation has a much larger amplitude of the semi-diurnal Kelvin wave than the free run and the Viking observations. This suggests that there is a resonant forcing associated with the analysis increments calculated with the 6-hr assimilation window. In order to test this hypothesis we recall that Wilson and Hamilton (1996) showed that zonal modulation of the sun-synchronous wave can be excited by (stationary) topographic forcing, as introduced in Section 2.5. An external forcing with wave number m could scatter the migrating tide (s=σ) into spatial harmonics σ−m and σ+m:
Equation (12) shows that m=4, corresponding to a wavenumber 4 topographic forcing, can indeed excite a wave of the form cos(−2λ+2t), or eastward-propagating semi-diurnal wave. Thus m=2 can yield a scattering of the diurnal westward (migrating) tides into A-1,1 (the diurnal east wave 1, or DE1). We are now interested in determining whether it is possible that the 6-hr assimilation cycle could introduce a geographically fixed forcing with a wave number 4 pattern that would also excite the non-migrating semi-diurnal east wave, or SE2. The equatorial temperature analysis increments in the original 6-hr TES assimilation do show a significant wave-4 pattern prevalent in several vertical levels, such as surface 800, 100–30 and 20–10 Pa (Fig. 2a). This introduces a forcing similar to a wave 4 topographic feature that excites SE2 [eq. (12)]. The use of 6-hr assimilation windows and a localisation length of 600 km produces a maximum amplitude in wave 4 (Fig. 2a), and thus maximum Kelvin wave forcing. This is because the analysis increments due to observations are computed four times a day, and the impact of the observations is constrained by the spacecraft orbit pattern and localisation length L=600 km, which keeps the analysis increments separated and limits their horizontal extent. The amplitude of the wave 4 in the analysis increments is significantly decreased by broadening the localisation to 1200 km (Fig. 2b and 3a). The wave 4 forcing in the analysis increments can also be modified by changing the updating hours, thus shifting the longitude phase of the pattern. Figure 2c shows the equatorial temperature analysis increment pattern after changing the assimilating hours to 03, 09, 15 and 21. It is clear that this wave 4 pattern shifts eastward compared to updating observation at hour 00, 06, 12 and 18, and the resulting resonance in Kelvin wave is also slightly diminished (Fig. 3b compared to Fig. 1d), in this case, the SE2 response is slightly diminished (Fig. 3b compared to Fig. 1d). Actually, the resulting SE2 will be a result of forcing by the ‘real’ physics and the spurious forcing. These two contributions may interfere constructively or destructively.
A 12-hr assimilation window with a short localisation of 600 km produces a wave 2 forcing (Fig. 2d), as could be expected from the previous arguments. Thus a 12-hr assimilation window induces a wave 2 pattern in the temperature analysis increment (Fig. 2d), serving as a spurious forcing to modify the wave number 1, eastward-propagating diurnal tide (DE1) (Fig. 4f). The combination of DE1 and the diurnal migrating tide yield an interference pattern having two maxima (Fig. 4c). It is obvious that the near-resonant Kelvin waves are sensitive to the DA implementation, while the migrating tides do not show this sensitivity. The non-uniqueness of the tide results show it is highly dependent on aspects of the DA system implementation (the assimilation phase and the localisation length scale), the first of which plays a more significant role by altering the forced and resonantly enhanced Kelvin wave considerably.
As noted previously in Section 2.5, the resonance condition is an intrinsic feature of Mars. In a resting atmosphere, the resonance will depend on atmospheric temperature, and the atmosphere is closest to resonance when the atmosphere is coldest, during the aphelion season (~Ls 100 in this study) (Wilson and Hamilton, 1996). Topography yields a way for forcing by the westward-propagating (migrating) tides to project onto eastward-propagating modes. The assimilation is able to create a spurious forcing that can be amplified by the resonance. Consequently, the amplitude and/or phase of the new response can be notably different from the ‘natural’ forcing, which is what we have seen, both for the 6-hr and the 12-hr updates. Changing the phase of the assimilation updates lead to changes in phase of the eastward-propagating components of the tide response, as we would expect. The migrating components are quite insensitive (figure not shown).
We have tested the hypothesis that large amplitude semi-diurnal and diurnal Kelvin waves in the analysis are spuriously forced by the use of 6-hr and 12-hr assimilation windows with short localisation scales (600 km), which concentrate the analysis increments into zonal wave number 4 and 2, respectively. In order to assess whether this spurious forcing can be eliminated (or at least reduced) we replaced the 6-hr assimilation window used in the LETKF with either a 2-hr or a 1-hr assimilation window. This reduction should actually improve the accuracy of the LETKF since shorter assimilation windows maintain a more linear growth of the forecast perturbations, as assumed in the development of the EnKF schemes (Kalnay et al., 2007). Compared to the 3D-LETKF (Fig. 5a, c and e), in which observations are gathered together in the centre of a 6-hr window, 4D-LETKF, in which observations are assimilated at their correct time of occurrence within a 1-hr window, yields more realistic observation increments. Note that since the observations are binned at hourly intervals, only 3D LETKF is conducted for 1-hr assimilation. This is particularly important for the Mars analysis, since the Martian boundary layer is dominated by large diurnal temperature variations and thus strong thermal tides (Greybush et al., 2012). For a forecast starting at hour 0, and using a traditional 6-hr 4D assimilation (Fig. 5b), a 9-hr forecast is launched first; observations are then combined together within 1-hr interval: t=3–3.5 hr, t=3.5–4.5 hr, t=4.5–5.5 hr, … and then t=8.5–9 hr. These observation files are then compared to forecasts at t=3, 4, 5, 6, 7, 8 and 9 hr to create an analysis at t=6 hr. Figure 5d and e show the basic structures of the 2-hr 4D and 1-hr 3D assimilation configurations used in this study. For the 2-hr 4D assimilation (Fig. 5c), also for a forecast starting at t=0 hr, observations are combined into 1-hr interval, i.e. t=1–1.5 hr, t=1.5–2.5 hr and t=2.5–3 hr. These observations are compared to forecasts at t=1, 2 and 3 hr, providing an analysis at t=2 hr (each model run produces 3-hr forecasts). For the 1-hr 3D assimilation (Fig. 5e) starting at t=0 hr, observations within t=0.5–1.5 are simply combined together to be compared to the forecast at t=1 hr, yielding an analysis at t=1 hr (each model run only generates a 1-hr forecast). The forecast RMSD difference between 4D and 3D LETKF are discussed in Section 4.5. For simplicity, in the following parts in this study, ‘6-hr/2-hr assimilation’ refers to ‘6-hr/2-hr 4D assimilation’ and ‘1-hr assimilation’ refers to ‘1-hr 3D assimilation’. Note that all runs of the MGCM (free runs, runs during assimilation and forecasts from analyses) in this study start and stop every hour to facilitate comparisons with observations.
Figure 6 reveals the tide amplitudes for the 2-hr and 1-hr TES assimilation results, both using the short localisation of 600 km. We can easily see that the amplitudes of A-2,2, or SE2, are significantly reduced compared to that of the 6-hr assimilation (Fig. 1d), and they even drop to a value below the free run (Fig. 1c). This is further demonstrated in the analysis increments, both in surface pressure (Fig. 7a–c) and 50 Pa temperature (Fig. 7d–f): weak wave 4 patterns exist in both the surface pressure analysis increment (Fig. 7a) and the 50 Pa temperature analysis increment (Fig. 7d) in the 6-hr TES assimilation, but not for the 2-hr (Fig. 7b and e) and 1-hr (Fig. 7c and f) TES assimilation. As a result, we confirm that the 6-hr (4 times/sol) assimilation leads to an erroneous geographically fixed forcing that excites the near-resonant semi-diurnal Kelvin wave to which the Mars atmosphere is particularly sensitive. This problem does not exist in short window assimilations like 2-hr and 1-hr windows, because the artificial separation of the analysis increments into four equally spaced groups disappears with short windows. In the following sections, we will evaluate how well the short window length analyses and their forecasts depict the state of the Martian atmosphere.
Observations at the two Viking Lander sites and the Curiosity site have provided surface pressure time series that could aid in the study of the seasonal variation of Martian climate. The Viking record consists of almost 4 Martian years, while Curiosity just completed a Martian year on June 24, 2014. Although the pressure observations have provided evidence for the Martian thermal tide full seasonal variation (Wilson and Hamilton, 1996; Bridger and Murphy, 1998), the current study is only focused on the aphelion season (~Ls 100).
Wilson and Hamilton (1996) and Bridger and Murphy (1998) conducted a harmonic analysis of the surface pressure tidal components at the two Viking Lander sites. Similar analysis was also conducted in this study. At Ls ~100, the corresponding Diurnal (S1) and Semi-diurnal (S2) amplitudes at the Viking Lander 1 (VL1) and Curiosity sites, as well as in all the experimental runs, are recorded in Table 1. Only modelling results are shown at the Viking Lander 2 (VL2), due to the lack of data during this particular season. As mentioned in Section 4.1, the surface pressure amplitudes at both VL1 and Curiosity site have also been normalised to compensate for atmospheric mass changes due to the CO2 condensation cycle and are shown in percentages of the local diurnal average surface pressure.
As shown in Table 1, at VL1, S1 varies little among different assimilation runs from 1.3 to 1.6%, and S1 values in all the experimental runs are stronger than in the observation (~1%). As VL1 is located at the fringe of the wave 4 interference pattern of S2 as shown in Fig. 1, S2 at VL1 is very sensitive to experiment details; the amplitude varies from 0.05% in the free run to 2.5% in the 6-hr TES assimilation, which is demonstrated to be unrealistic. S2 amplitude in all assimilation runs shows a much too strong signal at the VL1 location compared to in situ observations (0.3–0.4%). As suggested in Section 4.1, the S2 tide variation at the VL1 site depends on both A2,2 and A-2,2 and therefore is strongly related to the resonance in semi-diurnal Kelvin wave. Previous studies have found that under a uniform dust scenario, a strong Kelvin wave response during NH summer solstice (Ls ~100°) will result in a stronger diurnal tide at the Viking Lander sites (Wilson and Hamilton, 1996; Bridger and Murphy, 1998). Under the ‘seasonal dust’ scenario in which dust mixing ratio varies with season and latitude (Montmessin et al., 2004), the strong diurnal tide is suppressed at both VL1 and VL2 (Table 1). Due to the lack of high temporal resolution data in NH summer solstice at VL2, a strict comparison should only be made among modelling results, although we should also expect the consistent pattern of interference between eastward- and westward-propagating tides at this location. Observations by VL2 do support a small (less than 0.5%) S1 and S2 amplitude in this season, which confirms a seasonal decline in S2 amplitude observed at the other locations. Modelling results for VL2 also suggest S2 in a freely running model has a much smaller amplitude than in assimilation runs and TES 6-hr assimilation still shows a much stronger signal than in short window assimilations, while S1 amplitude varies little among different model runs.
Like VL1, the Curiosity site is also located in the fringe of the resonant wave 4 pattern, and thus sensitive to variations in the interference between A2,2 and A-2,2. Curiosity and VL1 are almost exactly 180° out of phase and hence are influenced by DE1 and SE2 in a very similar fashion (Fig. 1; Table 1). The interference pattern between the Kelvin and migrating tides will act roughly the same at each lander. Therefore, the increase in SE2 (and DE1) during NH summer solstice season leads to a decline in observed S2 at each lander. Compared to the S2 amplitude at the Curiosity tide (0.5–0.6%), despite a much too strong ‘artificial’ resonance in 6-hr TES assimilation, S2 amplitude in short window assimilations is still larger (0.9%). Although previous studies have shown that the observed normalised S2 within Gale is a reasonably close approximation of S2 in the surrounding area away from the crater, normalised S1 within the crater is larger than for the immediate surrounding area (Tyler and Barnes, 2013). An S1 amplitude of 3.1–3.3% in the assimilations in the present study should be reasonable given that the observation is ~3.4% and this is already an enhancement over the surrounding region (Table 1).
Aerosol radiative forcing (dust and/or water ice cloud) has a strong influence on the depth-weighted 50–10 Pa semi-diurnal tide amplitude. There is evidence suggesting that clouds, not dust, provide the dominant radiative forcing in non-dusty seasons, such as the aphelion season studied in this paper (Wilson et al., 2008; Kleinböhl et al., 2013; Steele et al., 2014). Tropical clouds contribute to strong daytime heating by absorbing upwelling IR radiation from the relatively hot surface. While clouds reduce nighttime surface cooling by IR emission at cloud base, there is also significant adiabatic heating in the 100–10 Pa region by the tide-induced downward vertical velocity, resulting in a net heating effect in the tropics (Wilson et al., 2008; Madeleine et al., 2012; Wilson and Guzewich, 2014).
Figure 8 presents the semi-diurnal tide amplitude of the temperature at 50 Pa. As discussed in the previous section and shown in Fig. 2a, the presence of a spurious wave 4 in the analysis increments with a 6-hr assimilation window leads to the excitation of a mid-level temperature semi-diurnal Kelvin wave (Fig. 8b and f) and is removed by using short window assimilations (Fig. 8c, d, g and h). It can be seen that for free runs, in the presence of RAC, the westward-propagating wavenumber 2, or A2,2, shows a significant signal in tropical regions (Fig. 8e). However, in the absence of RAC, the free run (Fig. 8a) shows a cold bias within this region, suggesting that adding ice cloud produces a substantial difference in the semi-diurnal tide amplitude through the region where the RAC is present (100~10 Pa). However, in the TES LETKF assimilations, the differences in A2,2 between with (Fig. 8f–h) and without RAC (Fig. 8b–d) are much smaller, suggesting that LETKF data assimilation can compensate for the absence of ice cloud parameterisation. This occurs because the assimilated observations implicitly contain information about heating from water ice clouds. Therefore, the assimilation run does ‘correct’ the cold tropical temperature bias in zonal mean temperature, as observed in the simulations here and in the MACDA assimilation (Wilson et al., 2008; Montabone et al., 2014). Moreover, without the presence of RAC (Fig. 8a–d), A2,2 is reasonably enhanced in the 2-hr and 1-hr LETKF assimilating TES profiles, while the free run and 6-hr assimilation may underestimate them (Fig. 3 in Kleinböhl et al., 2013, indicates that at 50 Pa, A2,2 should be around 2 K in the equator). The deficiency of the amplitude of A2,2 is largely due to the absence of RAC parameterisation in the free run and 6-hr assimilation, since A2,2 in the free run and 6-hr assimilation are largely enhanced with RAC parameterisation (Fig. 8e and f). This may further indicate that besides removing spurious resonance, short assimilation windows further ‘push’ the analysis closer to the true state even without the presence of necessary RAC parameterisation. This is consistent with the previous discussion that shorter assimilation window lengths are favourable for EnKF schemes like LETKF.
Since the twice-daily observations are insufficient to provide constraints on the semi-diurnal tides, it may not be possible to state that the short window assimilations are closer to the true state, or determine whether 1-hr assimilation indeed produces closer results to the true state than 2-hr assimilation. In fact, the observations are currently not very constraining (Wilson and Guzewich, 2014), which makes the identification of the optimal assimilation window length harder. However, we did show that the spurious resonance, apparent both in surface pressure and mid-level temperatures induced by the 6-hr assimilation cycle, could be removed in the short window assimilations. Moreover, the results with both 2-hr and 1-hr assimilation windows are quite similar, suggesting that 2-hr windows may be short enough to ensure linear error growth.
We next examine the impact of assimilation window length on the skill of forecasts (compared to future observations) at various lead times. In order to estimate the effect of assimilation window length on the Martian atmosphere prediction, as well as study the impact of RAC on the prediction, we also focus on the NH summer period, near the aphelion (Ls ~100) as in our previous tidal analysis. During this time period (NH summer), the dominant variability is from the tides. There is relatively low dust loading during this season and model errors are dominated by the aerosol vertical distribution, in which water ice clouds likely play the dominant role (Wilson et al., 2008, 2014). Little to no travelling waves and dynamical instabilities are present during this time period (Greybush et al., 2013). Figure 9 shows forecast RMSD with respect to temperature observation of 6-hr (4D, 3D), 2-hr (4D, 3D) and 1-hr (3D) assimilation under different scenarios, with and without RAC parameterisation during NH summer. Comparing 4D and 3D LETKF in 6-hr and 2-hr assimilations, 4D LETKF has smaller or almost equal RMSD than 3D LETKF, but the difference is not statistically significant (at a 5% level). RAC demonstrates again its significance in the MGCM: adding RAC yields a substantial reduction in the RMSD in both free running models and LETKF assimilations (the RMSD difference with and without RAC is statistically significant at a 5% level), by contrast, for the LETKF assimilations, RAC only reduces the averaged RMSD, by 2–5%. This indicates that the assimilated TES temperature observations already contain information about heating from ice clouds. With the presence of the RAC, the free run RMSD with observations is also much closer to the assimilation, about 20% higher rather than 80% as obtained without RAC. This would further suggest the importance of the RAC parameterisation, as indicated by Wilson et al. (2008), Madeleine et al. (2012) and Steele et al. (2014) that during the aphelion season the absence of a realistic representation of the radiative effects of the low latitudes water ice clouds may cause a significant systematic error in the Martian atmospheric prediction. This is confirmed by the tendency to have increasing errors right after the missing observations, especially in the absence of RAC, because during those times the ‘training’ of the model stops and the model is ‘pushed back’ towards the free run.
Although 1-hr assimilation shows higher frequency fluctuations, its forecast RMSD with observations is the smallest. However, considering that 1-hr/2-hr assimilation RMSD is comparing observations with 1-hr/2-hr forecasts and 6-hr assimilation RMSD is comparing observations with 6-hr forecasts, longer forecast lead times nearly always result in reduced prediction skill. Six-hour assimilation windows indeed show larger RMSD than 1-hr and 2-hr assimilation: the RMSD differences between 6-hr and short window assimilations, both with and without RAC, are statistically significant. However, the RMSD differences between 2-hr and 1-hr assimilations are not statistically significant whether we use RAC or not.
In order to choose the most appropriate assimilation window, a comparison of prediction that requires the same forecast length is needed. Six- to 336-hr (14-d) forecasts based on hour 00 and hour 12 analyses during sol 393–407 (Ls 96.8°–103.4°) are conducted and compared to observations within ±3 hr (e.g. 6-hr forecast is compared to observations within hour 3–9). Since there are fluctuations in representativeness error depending on where the orbit passes, averaged RMSDs over this time period (both with and without RAC) are used and shown in Fig. 10. The number of forecasts available to be used in the statistics is 29 for 24-hr forecast, 27 for 48-hr forecast, 25 for 72-hr … (decreasing by 2 each day, until 3 for 336-hr). Statistical error bars are also shown for each experimental run, indicating the errors are similar ranges among different assimilations runs.
For the forecasts including RAC parameterisation, it is clear that although 6-hr assimilation benefits from more observations (at hour 00, 6-hr LETKF assimilates observations from hour 21 (previous day) to hour 03; 2-hr assimilation uses observations from hour 21 (previous day) to hour 01, and 1-hr assimilation utilises observations from hour 21 (previous day) to hour 0.5), the different window lengths have similar predictive skill prior to the ‘saturation’ (~2 d). After reaching ‘saturation’, the free run and the assimilation forecasts have RMSDs that are close to each other, which means they lost all ‘predictability’. Using different observation windows of ±0.5 and ±1 hr for verification gave similar results (not shown here).
For model runs without RAC parameterisation, the 6-hr assimilation performs slightly better for short-range forecast (1–3 d) and short assimilation windows yield the smallest RMSDs for mid-range forecast (7–9 d), although their differences are not significant at a 5% level. In the absence of the RAC, ‘saturation’ occurs at day 9, much delayed compared to the ~2 d saturation found with RAC. The low ‘predictability’ (short time until error saturation) in the RAC scenario is due to the near-perfect nature of the freely running model including RAC. In this season with no active dynamics, getting the radiative heating correctly is the key issue.
From the above discussion, we may conclude that the short assimilation window lengths do not play a significant role in improving the Martian atmospheric prediction, at least in the aphelion (NH summer) period with low dynamical (chaotic) perturbation growth rates, although they do remove the spurious resonance in the Kelvin wave component of the tides as well as improve the semi-diurnal migrating tides in mid-level temperature simulation. Water ice cloud heating is an important factor that could affect the prediction in NH summer when aerosol loading is dominated by water ice clouds. When looking at the Martian atmospheric prediction, time to error saturation (RMSD from forecast matching free run) is an important index determined by both the radiative time scale and the dynamical error doubling time (Rogberg et al., 2010; Greybush et al., 2013).
Since every model has its own deficiencies, including both ‘external’ errors (e.g. a subgrid-scale parameterisation could never be perfect since it is hard to parameterise physical processes happening within the model grid) and ‘internal’ errors (imperfect initial condition, uncertainties and model instabilities), bias always exists within models. Model bias itself accounts for much of the temperature analysis error; for example, a Mars model requires parameterisation assumptions for water ice cloud and dust properties, and inaccuracies within these parameterisation processes will result in model bias from the true state. As a result, an empirical bias correction method based on Danforth et al. (2007) is applied to improve the model prediction. One of the advantages of data assimilation is to help compensate for imperfect model parameterisation by approximating the missing model forcings, treating the model as a black box and only considering forecast error statistics. These bias statistics can then be used to provide insights for improving model parameterisations and aerosol distributions, which are key goals for Mars science. The basic idea of empirical bias correction is to estimate model bias (in temperature, wind, surface pressure and surface stress) as the time-mean analysis increment over a period of time and then add the calculated bias to each of the ensemble member forecast before conducting the analysis. Greybush et al. (2012) conducted a sliding window bias correction centering at the analysis date using 10-sol time-averaged analysis increments from a previous assimilation run. In our experiments, we assume that there is no previously existing assimilation run, and the bias correction fields created from the analysis increments time averaged within the previous seven sols are used to correct the forecast at each time step – we call this as ‘online’ bias correction, which can save considerable computational time and resources. These results were compared to those using the bias correction method that relies upon a previous run, and their difference is only 1% in RMSD (figure not shown), which is small enough to justify using ‘online’ bias correction, since computation time is halved. Considering a 2-sol spin-up time from sol 388 (Ls 94.5°), we utilised the biases starting from sol 390 (Ls 95.4°) and started bias correction after day 391(Ls 95.9°). If the date that is 7 sols prior to the current analysis time is prior to sol 390, then the starting date to calculate the 7-sol average is hard-coded to sol 390, e.g. sol 391 uses daily-averaged bias of sol 390, and sol 392 uses daily-averaged bias averaged from sol 390 to 391. Since the model without RAC parameterisation is known to have large deficiencies in representing mid-level temperature, bias correction in the model that is far from perfect ought to have larger improvement.
Greybush et al. (2012) improved empirical bias correction by treating separate bias correction fields for different times of day (hour 00, 06, 12 and 18, for traditional 6-hr reanalysis), which is termed as ‘diurnal’ bias correction since different times of a day have different bias structure. Bias patterns at level 10 (near the top of TES observation levels) at different times of a day (hour 00, 06, 12 and 18) of 6-hr and 1-hr TES assimilation are shown in Fig. 11, indicating that bias structure varies at different times of a day. Daily average bias structures were shown in Fig. 7d and f; the difference between Figs. 7 and 11 illustrates the need for diurnally varying bias correction. Figure 11 also shows the corresponding observation path at different times of a day (hour 00, 06, 12 and 18) over 7 sols (sol 397–403, or Ls 98.6°–101.7°). It is also clear that within 6-hr window there are more observation tracks, and the observation path transition at the same hour (say hour 00) between different sols just shows a smooth position shift (Fig. 11a–d). However, since the 1-hr observation paths are more sporadic, the observations at the same hour (say hour 00) of different sols contain not only a significant position shifting, but also a significant contrast in observation numbers (Fig. 11e–h). Besides, the number of observations of different hours in the 1-hr observation path may vary much more than that in the 6-hr observation path. Due to the reasons presented above, larger sampling (representativeness) errors are expected in the 1-hr bias correction.
Table 2 clearly shows that in the absence of the RAC parameterisation, the ‘online’ sliding window diurnal bias correction runs have reduced the forecast RMSD compared to both no bias correction run and the run with ‘online’ sliding window daily average bias correction (difference between no bias correction run and daily average bias correction run, and between daily average bias correction run and diurnal bias correction run are both statistically significant at the 5% level). Note that the RMSDs of both the above-mentioned ‘online’ sliding window diurnal bias correction runs and ‘online’ sliding window daily average bias correction runs are lower than the runs without bias correction but with water ice clouds, passing the 5% significance test (Table 2; Fig. 9). We choose the forecasts at hour 00 and 12 during sol 397–407 (Ls 98.6°–103.4°) from both 1-hr and 6-hr ‘online’ sliding window diurnal bias correction runs, conduct the prediction runs with bias correction at their assimilation times and compare forecast (1–9 d forecast) with observations within ±3 hr, as described in Section 4.5. In this sense, we can interpret of the bias correction acting as ‘missing’ model physics during the forecast. The bias used in the prediction bias correction experiment is obtained from the sliding window 7-sol composite prior to the forecast time.
The bias-corrected 5-d (120 hr) observation increments suggest a big part of RMSD is persistent temperature bias in this season (Fig. 12b, c, f and g). For 1-hr bias correction, temperature differences are amplified in the upper levels near the SH polar front region (Fig. 12b and c), suggesting that the dynamics here interplay poorly with the 1-hr bias correction. This might be explained by the sensitivity of upper-level baroclinic waves to large sampling errors in 1-hr diurnal bias. The diurnal bias works better for both 1-hr and 6-hr corrections in the tropics (30°S–30°N) where baroclinic waves are not prevailing (Fig. 12b, c, f and g). Note that the observation increments in the bias-corrected runs in the absence of RAC (Fig. 12b, c, f and g) are also smaller than in the original runs without bias correction but with RAC parameterisation (Fig. 12d and h), suggesting that bias correction is indeed compensating for the absence of water ice clouds.
We examined the diurnal cycle of Mars atmospheric analyses created using the LETKF ensemble analysis system, the MGCM dynamical forecast model and the TES spacecraft observed temperature profiles. The Mars atmosphere is expected to be closest to resonance during the relatively cold aphelion (NH summer) season (Zurek, 1988; Wilson and Hamilton, 1996), therefore the Kelvin wave surface pressure response is most sensitive to forcing (realistic and/or spurious, e.g. from data assimilation) during this season. Kelvin wave resonance shows a large sensitivity to assimilation phase and moderate sensitivity to localisation radius, as discussed in Section 4.1. The use of a ‘traditional’ 6-hr assimilation window creates wave number 4 patterns in the daily averaged analysis increment, serving as an external forcing to trigger spurious resonance, which is manifested as a strong wave 4 pattern in surface pressure. The mechanism for generating this resonance is similar to how a wave 4 pattern in topography modulation excites an eastward-propagating semi-diurnal east wave. Harlim and Majda (2008) also noted a lack of observability when assimilating regularly-spaced, sparse observations with a resonant period in a simplified chaotic model; the data assimilation system cannot discern the true signal from among the various Fourier harmonics in the aliasing set based on observations alone.
In order to test whether this resonance mechanism is responsible for the large amplitude semi-diurnal tidal wave observed in our study, we recomputed the analysis using different assimilation window lengths. As expected, a 12-hr assimilation window led to a resonant wave 2 pattern in the analysis increments. Short assimilation windows such as 2 hr and 1 hr not only eliminated this artificial resonance induced by the traditional (for Earth systems) 6-hr assimilation, but also reduced the RMS difference between the forecast and the observations, and better compensated for the absence of the RAC parameterisation by considering the heating implied in the assimilated temperature profiles.
It is not surprising that the short assimilation windows show an advantage over the 6-hr window within LETKF. The ensemble Kalman filter assumes that the ensemble forecast perturbations with respect to the mean and the observational errors are both Gaussian, so that when they are combined, using the Kalman filter formulation, the analysis ensemble perturbations are also Gaussian (Hunt et al., 2007). This implies that if there are sufficient observations, windows that are short enough that the model forecast perturbations evolve linearly and thus remain Gaussian, are advantageous for EnKF. In our experiments with the assimilation of TES temperature profiles on the Mars atmosphere we have found that assimilation windows of 1 and 2 hr succeed in eliminating the resonance observed when using a 6-hr window, and their results are similar except that 1-hr analysis tends to have slightly larger (but not significantly different) RMS differences with the observations than 2-hr analysis in the absence of the RAC parameterisation. This suggests that for very short windows, since there are fewer observations per cycle, the analysis increments may be more sensitive to noise from observation error and sampling error in the background ensemble error covariance matrix, and that the optimal EnKF should be short enough to ensure linear growth, but not too short in order to ensure enough observations are processed simultaneously. By contrast, 4D-Var data assimilation is optimal for very long windows, in which the assimilation is exposed to the maximum possible number of noisy observations (Pires et al., 1996). They found, however, that for very long windows, due to the presence of multiple minima in the variational cost function, it was necessary to introduce a Quasi-static Variational Assimilation (QVA) in which the assimilation windows were started small, in order to reach the true single minimum, and then increased in small increments to reach the optimal long window. Singleton (2011) showed that for a toy coupled ocean-atmosphere model, the optimal coupled variational assimilation was reached using very long windows and the QVA algorithm. In contrast, the optimal coupled LETKF analyses were obtained with very short windows, but for optimal windows, the analyses errors were similar for the EnKF and the 4D-Var approaches.
Although all the assimilations showed significant improvements in forecast RMSD with respect to independent observations compared to a freely running model, short window assimilations further reduced the RMSD compared to the 6-hr assimilation. Considering that short window assimilations use short-term forecast to compare to nearby observations, this provides an advantage to short windows so that it may not be fair to compare among different assimilation forecasts with different assimilation window length. Therefore, we also launched 14-d forecasts during the NH summer (Ls~100), starting from hour 00 and 12 initial conditions, and calculated the averaged forecast RMSD compared to all the observations within ±0.5, 1 and 3 hr of the forecast time. With these longer forecasts, the assimilation window length appears to have only a minor impact on the forecast skill. Instead, the model forecast skill depends strongly on the role of water ice clouds (RAC) during NH summer when water ice cloud is prevalent. The ultimate goal of this study is to provide the basis for a multi-year, short assimilation window reanalysis that contains diurnal, seasonal and interannual variability, removes spurious resonance and gets as close as possible to the ‘true’ state of the Martian atmosphere.
To further improve model forecasts, we developed an ‘online’ sliding window bias correction method, which does not rely on a previous analyses. Since bias structure is different at different times of day due to the changing satellite track, ‘diurnal’ bias correction, which treats bias correction fields differently for different times of a day, was found to further reduce forecast RMSD as well as enhance prediction. It was shown that diurnal bias correction improves the prediction compared to conducting the daily average bias correction, but correcting the model bias every hour in the forecast induces larger sampling errors than correcting it using a longer window length (every 6 hr). In this sense, the bias error estimation in the mid-to-long range forecasts benefits from longer windows. The bias correction was found to be more effective than the water ice cloud in reducing the systematic model errors. Although the bias correction is obtained empirically and is not based on physical arguments, it can serve as guidance to modellers for the improvement of the physical parameterisations.
As discussed in Section 4.3, a significant issue is that observations twice per sol may not be sufficient to strongly constrain the tides in an assimilation. Since tides will be most dependent on the state of the aerosol field, one possible next step is to analyse temperatures and aerosols simultaneously, taking into account the correlations between aerosols and temperatures captured by the ensemble. Discussion of the impact of MCS profiles, including cross-track observations (with six local times), on features of a Mars reanalysis is part of our on-going work.
Our finding of assimilation-induced artefacts in atmospheric motions may have relevance to the Earth as well. Mitchell et al. (2014) highlights the dynamical similarities between the Martian atmosphere and the Earth's stratosphere. Tidal motions are important in the Earth's upper atmosphere as well (Hagan and Roble, 2001; Coy et al., 2011; Sakazaki et al., 2012). Swinbank et al. (1999) note a contrast in tidal amplitudes between simulations and assimilations for the Earth's stratosphere, suggesting possible improvement that should be made to the satellite retrievals. With the Earth's upper atmosphere being an increasing area of research for assimilation studies (e.g. Polavarapu et al., 2005; Sankey et al., 2007; Hoppel et al., 2013), the impact of satellite orbits and assimilation windows should be considered. However, the diurnal period is special to Mars and the comparable Kelvin resonance period on Earth is 33 hr, so the potential of forcing resonant atmospheric Kelvin modes is limited on Earth, as discussed in Section 2.5. Other normal modes might interact with the 6-hr assimilation interval in terrestrial data assimilation systems and this issue deserves detailed study.
31Note that 1 sol=1 Mars day=24 hr 37 min. A Martian hour=1/24 of a sol.
We acknowledge the NASA PDS for providing TES, Viking and Curiosity data. This research was supported by NASA grants NNX07AM97G and NNX11AL25G. We also thank John Harlim, François Forget, Thomas Navarro, anonymous reviewers, and members of the data assimilation and Mars atmosphere communities for insightful discussions that benefited this paper.
Christensen P. R. , Bandfield J. L. , Hamilton V. E. , Ruff S. W. , Kieffer H. H. , co-authors . Mars Global Surveyor Thermal Emission Spectrometer experiment: investigation description and surface science results . J. Geophys. Res . 2001 ; 106 : 23823 – 23871 .
Conrath B. J. , Pearl J. C. , Smith M. D. , Maguire W. C. , Christensen P. R. , co-authors . Mars Global Surveyor Thermal Emission Spectrometer (TES) observations: Atmospheric temperatures during aerobraking and science phasing . J. Geophys. Res . 2000 ; 105 : 9509 – 9519 .
Coy L. , Eckermann S. D. , Hoppel K. W. , Sassi F . Mesospheric precursors to the major stratospheric sudden warming of 2009: validation and dynamical attribution using a ground-to-edge-of-space data assimilation system . J. Adv. Model. Earth Syst . 2011 ; 3 : M10002 .
Eluszkiewicz J. , Moncet J. L. , Shephard M. W. , Cady-Pereira K. , Connor T. , Uymin G . Atmospheric and surface retrievals in the Mars polar regions from the Thermal Emission Spectrometer measurements . J. Geophys. Res . 2008 ; 113 : E10010 .
Forbes J. M . Tidal and planetary waves . The Upper Mesosphere and Lower Thermosphere: A Review of Experiment and Theory (eds. R. M. Johnson and T. L. Killeen) . 1995 ; Washington, DC : AGU . 67 – 87 .
Forbes J. M. , Hagan M. E . Diurnal Kelvin wave in the atmosphere of mars: Towards an understanding of “stationary” density structures observed by the MGS accelerometer . Geophys. Res. Lett . 2000 ; 27 : 3563 – 3566 .
Forget F. , Wanherdrick Y. , Lewis S. R . Validation of the Mars General Circulation Model and Climate Database with New Spacecraft Observations . 2001 . Work Package 7, Tech. Note 11369/95/NL/JG, European Space Agency, Paris .
Gómez-Elvira J. , Armiens C. , Castañer L. , Domínguez M. , Genzer M. , co-authors . REMS: The environmental sensor suite for the Mars Science Laboratory Rover . Space Sci. Rev . 2012 ; 170 : 583 – 640 .
Greybush S. J. , Kalnay E. , Hoffman M. J. , Wilson R. J . Identifying Martian atmospheric instabilities and their physical origins using bred vectors . Q. J. Roy. Meteorol. Soc . 2013 ; 139 : 639 – 653 .
Greybush S. J. , Wilson R. J. , Hoffman R. N. , Hoffman M. J. , Miyoshi T. , co-authors . Ensemble Kalman filter data assimilation of Thermal Emission Spectrometer temperature retrievals into a Mars GCM . J. Geophys. Res . 2012 ; 117 : E11008 .
Haberle R. M. , Gómez-Elvira J. , de la Torre Juárez M. , Harri A. M. , Hollingsworth J. L. , co-authors . Preliminary interpretation of the REMS pressure data from the first 100 sols of the MSL mission . J. Geophys. Res . 2014 ; 119 : 440 – 453 .
Hagan M. E. , Roble R. G . Modeling diurnal tidal variability with the National Center for Atmospheric Research thermosphere-ionosphere-mesosphere-electrodynamics general circulation model . J. Geophys. Res . 2001 ; 106 : 24869 – 24882 .
Hoffman M. J. , Greybush S. J. , Wilson R. J. , Gyarmati G. , Hoffman R. N. , co-authors . An ensemble Kalman filter data assimilation system for the martian atmosphere: implementation and simulation experiments . Icarus . 2010 ; 209 : 470 – 481 .
Hoppel K. W. , Eckermann S. D. , Coy L. , Nedoluha G. E. , Allen D. R. , co-authors . Evaluation of SSMIS upper atmosphere sounding channels for high-altitude data assimilation . Mon. Wea. Rev . 2013 ; 141 : 3314 – 3330 .
Lee C. , Lawson W. G. , Richardson M. I. , Anderson J. L. , Collins N. , co-authors . Demonstration of ensemble data assimilation for Mars using DART, MarsWRF, and radiance observations from MGS TES . J. Geophys. Res . 2011 ; 116 : E11011 .
Lee C. , Lawson W. G. , Richardson M. I. , Heavens N. G. , Kleinböhl A. , co-authors . Thermal tides in the Martian middle atmosphere as seen by the Mars Climate Sounder . J. Geophys. Res . 2009 ; 114 : E03005 .
Lewis S. R. , Read P. L. , Conrath B. J. , Pearl J. C. , Smith M. D . Assimilation of thermal emission spectrometer atmospheric data during the Mars Global Surveyor aerobraking period . Icarus . 2007 ; 192 : 327 – 347 .
Liu J. J. , Richardson M. I. , Wilson R. J . An assessment of the global, seasonal, and interannual spacecraft record of Martian climate in the thermal infrared . J. Geophys. Res . 2003 ; 108 : 5089 .
McCleese D. J. , Schofield J. T. , Taylor F. W. , Calcutt S. B. , Foote M. C. , co-authors . Mars Climate Sounder: an investigation of thermal and water vapor structure, dust and condensate distributions in the atmosphere, and energy balance of the polar regions . J. Geophys. Res . 2007 ; 112 : E05S06 .
Mitchell D. M. , Montabone L. , Thomson S. , Read P. L . Polar vortices on Earth and Mars: A comparative study of the climatology and variability from reanalyses . Q. J. Roy. Meteorol. Soc . 2014 ; 141 : 550 – 562 .
Montabone L. , Lewis S. R. , Read P. L . Interannual variability of Martian dust storms in assimilation of several years of Mars global surveyor observations . Adv. Space. Res . 2005 ; 36 : 2146 – 2155 .
Montabone L. , Lewis S. R. , Steele L. , Read P. L. , Ruan T. , co authors . Mars analysis correction data assimilation: a multi-annual reanalysis of atmospheric observations for the red planet . In: 4th World Climate Research Programme International Conference on Reanalyses . 2012 . Maryland, Silver Spring .
Montmessin F. , Forget F. , Rannou P. , Cabane M. , Haberle R. M . Origin and role of water ice clouds in the Martian water cycle as inferred from a general circulation model . J. Geophys. Res . 2004 ; 109 : 10004 .
Navarro T. , Forget F. , Millour E. , Greybush S. J . Detection of detached dust layers in the Martian atmosphere from their thermal signature using assimilation . Geophys. Res. Lett . 2014 ; 41 : 6620 – 6626 .
Rogberg P. , Read P. L. , Lewis S. R. , Montabone L . Assessing atmospheric predictability on Mars using numerical weather prediction and data assimilation . Q. J. Roy. Meteorol. Soc . 2010 ; 136 : 1614 – 1635 .
Sakazaki T. , Fujiwara M. , Zhang X. , Hagan M. E. , Forbes J. M . Diurnal tides from the troposphere to the lower mesosphere as deduced from TIMED/SABER satellite data and six global reanalysis data sets . J. Geophys. Res . 2012 ; 117 : D13108 .
Szunyogh I. , Kostelich E. J. , Gyarmati G. , Kalnay E. , Hunt B. R. , co-authors . A local ensemble transform Kalman filter data assimilation system for the NCEP global model . Tellus A . 2008 ; 60 : 113 – 130 .
Wilson R. J. , Millour E. , Navarro T. , Forget F. , Kahre M. A . GCM simulations of aphelion season tropical cloud and temperature structure . Mars Atmosphere: Modeling and Observations, 5th International Workshop . 2014 ; UK : Oxford .
Zhang K. Q. , Ingersoll A. P. , Kass D. M. , Pearl J. C. , Smith M. D. , co-authors . Assimilation of Mars Global Surveyor atmospheric temperature data into a general circulation model . J. Geophys. Res . 2001 ; 106 : 32863 – 32877 .