Assimilation of Both Column- and Layer-Integrated Dust Opacity Observations in the Martian Atmosphere.

A new dust data assimilation scheme has been developed for the UK version of the Laboratoire de Météorologie Dynamique Martian General Circulation Model. The Analysis Correction scheme (adapted from the UK Met Office) is applied with active dust lifting and transport to analyze measurements of temperature, and both column-integrated dust optical depth (CIDO), τ ref (rescaled to a reference level), and layer-integrated dust opacity (LIDO). The results are shown to converge to the assimilated observations, but assimilating either of the dust observation types separately does not produce the best analysis. The most effective dust assimilation is found to require both CIDO (from Mars Odyssey/THEMIS) and LIDO observations, especially for Mars Climate Sounder data that does not access levels close to the surface. The resulting full reanalysis improves the agreement with both in-sample assimilated CIDO and LIDO data and independent observations from outside the assimilated data set. It is thus able to capture previously elusive details of the dust vertical distribution, including elevated detached dust layers that have not been captured in previous reanalyzes. Verification of this reanalysis has been carried out under both clear and dusty atmospheric conditions during Mars Years 28 and 29, using both in-sample and out of sample observations from orbital remote sensing and contemporaneous surface measurements of dust opacity from the Spirit and Opportunity landers. The reanalysis was also compared with a recent version of the Mars Climate Database (MCD v5), demonstrating generally good agreement though with some systematic differences in both time mean fields and day-to-day variability.

Introduction

The dust cycle is a key component of the Martian climate, and is extremely important for understanding the interannual, seasonal and synoptic evolution of the Martian environment. (e.g., Kahre et al., 2017; Newman et al., 2002a, and references therein). Intensive measurements of atmospheric temperature and dust extending over more than ten Mars years (MY) now exist with unprecedented spatial coverage, thanks to various orbital spacecraft. Such observations have already helped to improve our understanding of Mars' weather and climate. However, the incomplete coverage of these measurements across the planet constrains our ability to study the general circulation in full detail, particularly those aspects related to dust opacity. For instance, the Thermal Emission Imaging System (THEMIS) carried by the Mars Odyssey (MO) spacecraft can provide multi‐annual measurements of Column Integrated Dust Opacity (CIDO), but its coverage in space and time is quite limited.

On the other hand, numerical models provide four‐dimensional simulated data with moderate to high temporal and spatial resolution and complete coverage in space and time, but often fail to reproduce the dust cycle's full range of variability. Various authors, starting with Newman et al. (2002b) (see also Kahre et al. (2017) for a review and Gebhardt et al. (2020b) and Gebhardt et al. (2020a) for more recent studies) showed that a global circulation model (GCM) could capture the onset and growth of regional dust events, but did not realistically capture the observed interannual variability. In particular, they could not reproduce the relatively “quiet” year of dust activity that occurs immediately after a simulated global dust storm (GDS) year. Others have sought to take additional factors into account, such as the finite extent of the surface dust reservoir (e.g., Pankine & Ingersoll, 2004; Szwast et al., 2006) or nonlinear effects associated with the “shadowing” of pockets of dust behind rocks and boulders (Mulholland et al., 2013). But even the most sophisticated free‐running GCMs still struggle to capture realistic interannual variability associated with dust lifting and transport.

To aid in this task, data assimilation has become an optimal approach to provide a solution that is consistent with both observations and modeled physical constraints. Data assimilation corrects model‐predicted variables toward observations such that the resulting solution can represent the full observed variability of the climate. Such an assimilated record is often termed a “reanalysis” by analogy with the practice in Earth weather and climate forecasting, in which routine meteorological observations, collected primarily in order to initialize numerical weather forecasts, are “reanalyzed” sometime afterward over long periods using a uniformly consistent model and assimilation scheme to produce a self‐consistent climate record. This approach has been widely used as an effective tool in operational weather forecasting systems or climate models for analyzing meteorological variables for the Earth (e.g., Kalnay, 2003; Lorenc et al., 1991). This approach has already been used for a number of years to investigate tracer/chemical evolution in the Earth's atmosphere (Benedetti et al., 2018; Collins et al., 2001; Schutgens et al., 2010; J. Wang et al., 2004). Collins et al. (2001), for example, used an optimal interpolation approach to assimilate satellite retrievals of total column aerosol optical depth (AOD) over the Indian Ocean, which reproduced the daily variations of AOD at a single model grid point. J. Wang et al. (2004) used nudging to assimilate AOD into a nonhydrostatic atmospheric model, which captured the observed evolution of a dust event near Puerto Rico. More recently, Schutgens et al. (2010) applied the Local Ensemble Transform Kalman filter to assimilate AOD from the AERONET global surface observation network, which captured the evolution of AOD and also reduced uncertainties in model estimates of the evolving aerosol distribution. At the time of writing, around five major operational centers around the world use a variety of assimilation techniques, including optimal analysis (similar to the Analysis Correction scheme presented here), variational methods or ensemble Kalman filters to analyze observations of dust and aerosols from various sources (e.g., Benedetti et al., 2018, and references therein), such as AOD derived from orbiting or surface‐based platforms.

Few centers have assimilated dust profile observations, however, preferring instead to focus on achieving high horizontal resolution utilizing the much more abundant AOD measurements. Limited publications to date include the work of Yumimoto et al. (2008), who assimilated vertical profiles of the dust extinction coefficients in a regional dust transport model. In their study, the data from a ground‐based lidar network were interpolated to the vertical model levels for analyzing the model prognostic dust variables. More recently, Sekiyama et al. (2010) directly assimilated the total attenuated backscattering coefficient from the Cloud‐Aerosol Lidar and Infrared Pathfinder Satellite Observations (CALIPSO) mission into a global chemistry‐transport model. The measurements were averaged approximately to the model horizontal and vertical resolution before the assimilation. This approach has recently been extended by Cheng et al. (2019), who also assimilated CALIPSO profiles of aerosol optical depth using a 4D Ensemble Kalman Filter approach.

Data assimilation has also been applied to the Martian atmosphere with success. Lewis and Read (1995) implemented the analysis correction (AC) scheme (Lorenc et al., 1991) in a simple version of a Mars General Circulation Model (MGCM) in order to assimilate temperature profiles from the Pressure Modulator Infrared Radiometer (PMIRR) instrument on‐board the short‐lived Mars Observer spacecraft (1993). Their results showed that assimilation of such observations was feasible and that it improved the agreement between model and observations. Lewis et al. (2007) extended this approach to include dust tracer assimilation, which was combined with a full MGCM to assimilate thermal profiles and CIDO rescaled to a reference level (hereafter τ ref) using retrievals from the Thermal Emission Spectrometer (TES) on‐board MGS (M. D. Smith et al., 2003). The performance of the data assimilation system was validated against independent radio occultation measurements by Montabone, Lewis, Read, and Hinson (2006). This showed that combined temperature and τ ref assimilation was able to reduce discrepancies between the model and radio occultation data below 20 km, especially when dust amounts were large and changing rapidly, although some large discrepancies remained due to known inconsistencies between TES temperature profiles and radio occultation data. This approach was further extended by Holmes et al. (2020) to include assimilation of column dust optical depth measurements derived from Mars Reconnaissance Orbiter (MRO) MCS retrievals, together with measurements of water ice and ozone, into a version of the LMD/UKMGCM that also advects dust and other tracers with the analyzed winds.

An alternative approach to assimilation of Mars observations was developed by Hoffman et al. (2010) based on a complementary method using the ensemble Kalman filter (EnKF, Evensen, 2003) to assimilate TES temperature retrievals. They found generally improved agreement with TES temperature observations over a free‐running model, and in Greybush et al. (2012) the joint assimilation of TES temperatures with forcing using a 2D CIDO dust field from TES was shown to improve the agreement of model and TES temperatures further. This method has also been extended to include assimilation of column dust optical depths from both MGS/TES and MRO/MCS observations (Greybush et al., 2019) and a data set is publicly available known as EMARS.

The approach used by Lewis et al. (2007) first assimilated TES temperature profiles and τ ref without explicitly advecting the dust tracer field, using an empirical relation (Conrath, 1975) to prescribe the vertical distribution of dust. This system has subsequently been applied in several studies of Martian weather and climate (Lewis & Barker, 2005; Lewis et al., 2007; 2016; Montabone et al., 2005; Montabone, Lewis, Read, & Withers, 2006; Rogberg et al., 2010; Wilson et al., 2008), and both a three‐year reanalysis covering MY 24–27 using this assimilation system and another covering much of the MCS period (MY 28–32) have been published (Holmes et al., 2020; Montabone et al., 2014). Navarro et al. (2014) also assimilated Mars Climate Sounder (MCS) temperature profiles and modified the dust vertical distribution using its correlation with temperature. However, it is essentially different from the work presented here, in which dust observations are directly assimilated. The more recent study of data assimilation issues on Mars by Navarro et al. (2017) is more similar to the present work in including some cases that assimilated MCS dust opacity profiles using the EnKF method. Their study indicated some promise for this approach, although they only analyzed a part of MY 29 and noted some difficulties in capturing the diurnal variation in dust vertical distributions.

The data assimilation system developed here is based on the scheme described by Lewis et al. (2007). However, that scheme does not assimilate a vertically resolved dust distribution, only TES nadir retrievals of CIDO, and the model does not transport the dust actively. The newly available data set from MCS on board MRO (Kleinböhl et al., 2009) does provide vertically resolved, global measurements of the atmospheric dust distribution. With this new data set in hand, here we update the existing data assimilation system to better represent the Martian dust cycle. In later work we will use this to study the formation and life cycles of regional and global dust storms in detail.

Section 2 describes the Mars GCM and current data assimilation scheme, and the observations of Martian dust are described in Section 3. We outline how the assimilation was adapted and extended to incorporate dust profile observations alongside column dust opacities in Section 4. Sections 5 and 6 present verifications against in‐sample and out‐of‐sample observations respectively, while Section 7 describes a systematic comparison of the Mars Climate Database against the reanalysis. We conclude in Section 8.

Overview of Mars GCM and Data Assimilation Scheme

In this work the model used is based on the UK version of a three‐dimensional Martian Global Climate Model (UK‐LMD MGCM, v5.1.3) (Forget et al., 1999; Mulholland et al., 2013). The model combines a spectral dynamical solver at triangular truncation T31, corresponding to a 96 × 48 longitude‐latitude grid in real space, a tracer transport scheme and dust lifting and deposition routines, along with a full range of physical parameterizations.

The equations for a hydrostatic, adiabatic and inviscid gas surrounding a rotating spherical planet are cast in vorticity‐divergence form. In the vertical, levels are defined in terms of the terrain‐following σ coordinate system using a standard finite difference approach. There are 25 levels with the first three at 4, 19, and 44 m above the surface, to resolve detailed surface processes represented in the model. The model top varies in altitude over time but is typically at around 100 km, with a sponge layer (applying a linear drag on eddy vorticity and divergence) in the uppermost three levels to reduce spurious reflections of vertically propagating waves. There are typically 480 dynamical and 96 physics timesteps per sol (where a sol is a mean solar day on Mars).

The radiative transfer scheme calculates atmospheric absorption and emission due to carbon dioxide and airborne dust; the radiative effects of water vapor and ice are not included since our focus here is on the dust cycle. We rely, therefore, on the temperature assimilation to account for the radiative effects of clouds. The balance between incoming radiative flux and thermal conduction in the soil contributes to changes in surface temperature, using a surface thermal inertia field derived from TES and Viking (Forget et al., 1999) and topography from the Mars Orbiter Laser Altimeter on MGS (D. E. Smith et al., 2001). The surface roughness length z 0 is based on a global map compiled by Hébrard et al. (2012), and implemented in the UK‐MGCM by Mulholland et al. (2015).

The dust transport scheme is fully described by Newman et al. (2002a) and Mulholland et al. (2013) and includes dust lifting parameterizations (using constant wind stress thresholds), tracer advection, gravitational sedimentation and dry deposition. We assume a 1.5 μm particle size for simplicity, based on Mars Exploration Rover (MER) observations (Lemmon et al., 2004). The two most important distinct mechanisms responsible for the injection of dust into the atmosphere are thought to be dust lifting by near‐surface wind stress, and dust lifting by dust devils (Newman et al., 2002a).

The data assimilation scheme is based on the analysis correction sequential estimation (AC) scheme (Lorenc et al., 1991) but with modifications specific to Mars (Lewis et al., 2007). The assimilation step is computationally inexpensive compared with the rest of the model, and so is usually performed at each dynamical timestep (typically of 3 min). Lewis et al. (2007) describe the scheme in full detail. Temperature assimilations are the same as in that work, except for the observational data set used. They assimilated dust CIDO observations without advecting the dust tracer, instead setting the vertical distribution of dust opacity using an empirical relation following Conrath (1975). In this work we extend the dust assimilation to incorporate advective transport of radiatively active dust in the simulation model as well as to assimilate both CIDO and LIDO (Layer integrated dust opacity) observations; this is described in Section 4.

The ratio of observational error to first guess error used in the normalization factor (Lorenc et al., 1991, Equation 3.20) is set to 1 for assimilation of MCS temperature observations (as previously done for TES by Lewis et al., 2007), implying that the background forecast and observation errors are comparable. Following the study of ice opacity assimilation by Steele et al. (2014), we also set this ratio to 1 for the dust observations.

Observations of Martian Dust

Thanks to various spacecraft in orbit around Mars since 1997, measurements of atmospheric temperature and dust exist covering the Martian atmosphere over more than ten MYs. The instruments on board these spacecraft for determining temperature and dust in the Martian atmosphere that have been used for assimilation, such as the present study, include (amongst others) TES on MGS (M. D. Smith, 2004), THEMIS on MO (M. D. Smith, 2009) and MCS on MRO (Kleinböhl et al., 2009).

TES and THEMIS dust retrievals contain CIDO data only, while MCS data contain satellite observations (MCS v3 was used here for this initial proof of concept) with vertically resolved, asynoptically sampled global retrievals of atmospheric profiles of temperature and LIDO (McCleese et al., 2010). This contains information on the day‐to‐day variability of Martian weather from the near surface to the top of the middle atmosphere around 80 km altitude (Kleinböhl et al., 2009). The spacecraft have different operational periods and orbits, so their retrievals have different temporal and spatial coverage. Only THEMIS has overlapping observational periods with the other two datasets. Further details of the THEMIS data set, including the retrieval algorithm, can be found in M. D. Smith et al. (2000, 2003) and M. D. Smith (2009).

The analysis in this paper focuses on part of the MCS mapping period from MY28 L  = 110° (solar longitude) to MY29 L  = 330°. During this period, THEMIS CIDO data and MCS LIDO data are both available, although the MCS data during the period in MY29 between L  = 180° and L  = 270° were taken in limb staring mode and are considered by the MCS team to be poorly calibrated. In practice this means that many profiles in that period were missing data from deeper levels though still with reasonable horizontal coverage. The potential impact of this is considered in the out of sample validation below. The spatial coverage of THEMIS retrievals varies significantly with L , while MCS retrievals are more consistent and uniform except during the global‐scale dust storm (GDS) season (roughly from MY28 L  = 270°–305°; see Figure 1b). Coverage was restricted to the very early stage of the storm and poleward of 45°N throughout. The spatial coverage of these two datasets within the study period is shown in Figure 1.

Figure 1

Spatial and temporal distribution of available dust opacity data from Thermal Emission Imaging System and Mars Climate Sounder during the study period. The color scales show the number of measurements in 5° L and 3° latitude bins.

Because the dust is assumed by the THEMIS dust opacity retrieval algorithm to be well mixed, the THEMIS CIDO data is commonly reported rescaled (i.e., as τ ref) to a reference pressure of 610 Pa (M. D. Smith, 2009), to remove the effects of variable topography. However, this assumption may introduce uncertainty when the dust is not well mixed. When an intense detached dust layer exists (Heavens et al., 2011), for example, rescaling under the well‐mixed assumption could lead to an overestimate in τ ref.

It is also important to note that THEMIS dust observations are provided as an infrared absorption optical depth, while the modeled τ ref is the visible extinction optical depth. Using numerical experiments, M. D. Smith (2009) determined the conversion between IR absorption and visible extinction optical depth to be γ ∼ 1.3, and Clancy et al. (2003) used a scaling factor ɛ = 2 to convert from IR to visible for dust of size 1.5–2.0 μm. Lemmon et al. (2004) compared visible optical depths with MER measurements at 9 μm wavelength, and found agreement with Clancy et al. (2003). For simplicity, in this work the dust particle size is approximated by a constant 1.5 μm, which is reasonably consistent with various observational studies (Clancy et al., 2003; Lemmon et al., 2004; Pollack et al., 1995). Hence the conversion factor from THEMIS IR absorption optical depth to a model‐compatible visible extinction optical depth is 2.6, and in this work we imply the visible extinction optical depth when referring to CIDO and τ ref, unless otherwise stated. For MCS LIDO retrievals, the infrared opacities (at a wavelength centered on 21.6 μm) were multiplied by a factor of 7.3 to convert them to a visible equivalent (Montabone et al., 2015).

For a fully independent validation of the analysis, upward‐looking surface observations provide a bottom‐up view of CIDO that is independent of satellite‐based datasets, although only over particular locations. The MER missions, Spirit (14.57 °S, 175.48 °E) and Opportunity (1.95 °S, 5.53 °W) provide almost continuous data coverage during MY28 and MY29, concurrent with the study period in this work. Both rovers carried a Pancam camera, which included solar filters at 440 and 880 nm wavelengths. The effective dust particle radii based on the rovers' observations were 1.47 ± 0.21 μm for Spirit and 1.52 ± 0.18 μm for Opportunity (Lemmon et al., 2004). However, the 440 nm filter is significantly affected by a red leak (Lemmon et al., 2015). As there should not be a significant difference between the measured CIDO at 880 nm and CIDO at 700 nm (used as the model visible wavelength), we therefore rescale the CIDO measurements at 880 nm to the reference pressure 610 Pa, and make the comparison directly with the modeled sol‐averaged τ ref (also rescaled to 610 Pa).

About 10% of the THEMIS and MCS data are excluded from the assimilation and used for out‐of‐sample validation. The withheld data were taken as every tenth THEMIS data point and every tenth MCS vertical profile. Withholding this small fraction of the data set should not greatly affect the assimilated results. One should not be surprised that these validation data may be correlated with the assimilated data, which does weaken their usefulness in validating the model to some extent. However, these datasets were the best available for the assimilation itself.

Dust Data Assimilation With Active Transport

Initial conditions for the prognostic variables and dust tracers were taken from a free‐running spin‐up run, which was run for two years prior to the start of the assimilation. Temperature data from MCS was included in each assimilation (though not the free running model during spin‐up) using the method described by Lewis et al. (2007), and such that temperatures were always assimilated before the dust assimilation. An identical free‐running simulation without any assimilation, but with a fully active dust lifting and transport cycle tuned to reproduce plausible seasonal variations of dust loading (cf. Mulholland et al., 2013; Newman et al., 2002a), was run in parallel.

Previous assimilation studies using the UK‐LMD MGCM excluded active dust transport, instead just correcting the temperature profiles and τ ref. The dust distribution remained static in the absence of dust observations and the vertical distribution was prescribed using the Conrath (1975) distribution. In the new scheme presented here, the data assimilation system is updated to include full dust transport, lifting and sedimentation while correcting τ ref and/or the model's vertical dust distribution with observations.

CIDO Assimilation Only

In this configuration only the CIDO retrievals are assimilated. The sequence of operations is shown in Figure 2. First the dynamics timestep is integrated, and then the dust is advected to obtain the three dimensional (3D) distribution of dust mass mixing ratio q(x, m). The CIDO at position x, predicted from this distribution (τ C) is obtained by linearly summing up the layer‐integrated dust opacity (LIDO) within each model layer (τ LIDO). In the model the LIDO at horizontal position x in model level m is given by where q ext = (3Q ext)/(4ρr), q(x, m) is the dust mass mixing ratio, g = 3.72 m s−2 is the gravitational acceleration, Q ext is the extinction coefficient, ρ = 2500 kg m−3 and r = 1.5 μm are the density and radius of dust particles, respectively. p (x) is the surface pressure, and Δσ(m) is the layer thickness. The reference dust opacity τ ref at a reference pressure p ref is then determined from the model by

Figure 2

Sequence of operations in the new data assimilation scheme with active dust transport. Green boxes show initial conditions and input observations, blue boxes show data generated by Mars General Circulation Model (MGCM) integration, and red boxes show individual MGCM modules. Text in black applies to both column integrated dust opacity (CIDO) and layer‐integrated dust opacity (LIDO) assimilation, text in blue applies to CIDO assimilation only, and text in red applies to LIDO assimilation only. Only variables related to the data assimilation scheme are included.

The reference pressure p ref is arbitrary; to compare the results with observations, modeled CIDO values are rescaled to 610 Pa.

The advected dust opacity and mass mixing ratio fields are then used to integrate the physical parametrizations. Finally, τ ref and T are updated using data assimilated by the AC scheme, followed by increments to u and v in thermal wind balance. The dust transport scheme transports a 3D dust mass mixing ratio field, but the assimilated CIDO dust observations constrain τ ref only, so the dust mass mixing ratio at each model layer must be adjusted after the assimilation. This adjustment simply consists of a multiplicative scale factor λ(x), which ensures that the shape of the vertical dust profile at each horizontal grid point remains the same before and after the assimilation of τ ref: where variables without and with primes are before and after assimilation, respectively. Since the extinction coefficient and layer thickness are constant within a particular time step, Equation 1 implies that the dust mass mixing ratio q(x, m) is proportional to τ LIDO(x, m), and therefore the adjustment λ(x) can be applied directly to q(x, m). A similar assumption was also used when assimilating AOD on Earth (e.g., Collins et al., 2001; J. Wang et al., 2004).

LIDO Only

A more advanced method is required to make proper use of the MCS vertically resolved dust profiles. This section describes how retrievals of dust profiles from MCS are assimilated into the model by themselves (i.e., without assimilating CIDO). Figure 2 shows the procedure for assimilation of LIDO, and it is very similar to CIDO‐only, but now LIDO (τ LIDO) is analyzed directly. When the model layers have a smaller vertical spacing than the MCS measurements (typical in the lower and middle atmosphere), this approach avoids the direct interpolation of observational data to the model levels.

Since MCS does not (in general) take data at the same levels as those used in the model, we need to pre‐process the observed dust distribution. Our approach differs in approach from both Yumimoto et al. (2008) and Sekiyama et al. (2010), who both used the CALIPSO satellite LIDAR measurements of dust opacity in the Earth's atmosphere. Here, MCS dust retrievals are reported as dust opacities at atmospheric pressures typically 1–1.5 km apart, but their intrinsic vertical resolution is about 5 km (Kleinböhl et al., 2009), so the data set oversamples the actual MCS measurements. First we integrate the MCS dust retrievals vertically with a 5 km grid spacing in order to recover the observed LIDO. This ensures that the assimilated data has the same vertical resolution as the actual measurements. This preserves smaller‐scale vertical variability in the modeled dust profile that would be unresolved in the MCS observations.

The approach used here resembles the assimilation of thermal profiles into the UK‐LMD MGCM (Lewis et al., 2007). First, we use the modeled dust opacities τ LIDO(m) to predict the dust opacity within observation layer i. Model layers that overlap more than one observed layer are split linearly in ln p among the observed layers. Figure 3 shows an example where three model layers overlap one observed layer. In this instance the modeled dust opacity of the observed layer is

Figure 3

Schematic showing the calculation of the dust opacity increment Δτ LIDO(i) within observation layer i, given the observed layer‐integrated dust opacity within that layer τ LIDO(i) (right) and the modeled dust opacities τ LIDO(m) at the overlapping model levels (left).

The LIDO increment within this observation layer is then

From this, the LIDO increment at each model layer m due to observation layer i is

The horizontal assimilation then uses these increments, summed over contributions from each observation level i, to update the modeled dust field, following the standard procedure for τ ref.

Dust is transported in terms of dust mass mixing ratio, so the assimilation needs to correct this quantity. As dust mass mixing ratio is proportional to LIDO, it is multiplied by a factor , where the primed and unprimed quantities are the corrected and uncorrected LIDO values.

Joint CIDO and LIDO

To take advantage of both datasets simultaneously, we can assimilate both CIDO and LIDO together. In principle, one could use the measured dust profiles to correct the dust in model layers where there are observations, then use the CIDO data to correct the rest of the column. This avoids unnecessarily adjusting the vertical distribution using CIDO when part of the distribution has already been corrected using LIDO data. However, THEMIS and MCS measurements are not normally taken at the same time and place, so it is difficult to use both simultaneously at one location. Therefore we instead assimilate the LIDO and CIDO datasets independently.

Verification I.: In‐Sample Observations and Free‐Running Model

The methods described above were used to analyze various combinations of THEMIS and MCS observations obtained during Mars Years 28 and 29, representing a typical pair of years that include dusty seasons both with and without a planet encircling event. In this section we present results that compare assimilated analyses with a free‐running model simulation with full dust transport and seasonal variability and evaluate the convergence of the assimilation toward the input data. Further results and figures can be found in Section S1 of the Supporting Information S1.

The model‐predicted and assimilated values of τ ref (rescaled to 610 Pa) from each variant of the scheme were interpolated to the positions of THEMIS CIDO measurements. Figure 4 (red line) shows the global mean τ ref of these interpolated data over the course of the study period. This shows that all three of the reanalyzes converge to the assimilated THEMIS data outside the GDS period in MY28, but in contrast to the other variants, the LIDO‐only assimilation overestimates the peak τ ref during the GDS and misrepresents the timing of its onset. The free‐running simulation, on the other hand, captures some of the variability, but completely misses the development of the GDS around MY28 L  = 300°.

Figure 4

Global mean τ ref (rescaled to 610 Pa) over the study period for the free‐running simulation (green), the column integrated dust opacity (CIDO)‐only reanalysis (red), layer‐integrated dust opacity (LIDO)‐only reanalysis (cyan), and joint CIDO/LIDO reanalysis (magenta). THEMIS CIDO observations used in the assimilation are shown as triangles. Each point is an average over five sols.

Assimilating CIDO Only Versus MCS Dust Observations

It is also useful to compare the dust reanalysis with the observed time‐zonal mean dust distribution. A set of vertical dust distributions (designated MCS‐binned observations hereafter) were produced by sampling MCS dust profiles in 5° horizontal grids during daytime (local time 06:00–18:00) and nighttime (local time 18:00–06:00), after binning the data in L  = 5° intervals. The model results were interpolated to the same grid and averaged over the same L time windows, and restricted to altitudes where MCS‐binned observations were available before taking the zonal mean.

A comparison is shown in Figure 5, which shows results for two cases (a) at L  = 352.5° of MY28, close to the northern Spring equinox and (b) at L  = 122.5° in late northern Summer, also of MY28. Note the lack of observations in Figure 5a near the northern winter pole, which is due largely to an inability of MCS retrievals to distinguish suspended dust from CO2 cloud aerosol, so observations are largely filtered out from the analysis. With CIDO assimilation, the top of the dust layer when detached dust layers are absent is broadly similar to observations in the case shown, as it is in the free‐running model (Figure 5a middle frame). This likely indicates that the dust distribution at this time happened to be indistinguishable from the dust climatology in the model. However, elevated detached dust layers, such as seen in Figure 5b, cannot be reproduced in either the CIDO‐only reanalysis or the free‐running model (or any other model without an explicit parameterization of “rocket” dust storms (cf. C. Wang et al., 2018)). Such detached dust layers were observed in MCS night‐time retrievals (Heavens et al., 2011) and later confirmed by other instruments (Guzewich et al., 2013; M. D. Smith et al., 2013). In the version of the model used in this study, dust tends to be lifted to a lower height than the observed detached dust layer, and it is then well mixed all the way to the ground. Hence a successful reanalysis is likely to require assimilating vertically resolved dust measurements (i.e., LIDO) to reproduce the detached dust layers in a reanalysis.

Figure 5

Night‐time (18:00‐06:00 local time) zonal‐time mean dust opacity (km−1) during MY28 (a) without and (b) with detached dust layers. Time averages are over 5°L .

Assimilating LIDO Only Versus THEMIS Dust Observations

In the LIDO‐only assimilation, the period with a GDS is captured (Figure 4, cyan line) despite the limited MCS coverage during that period (Figure 1b). These observations are sufficient for the reanalysis to capture the initial condition and northern boundary condition of the GDS, but the inferred peak in τ ref is higher and later than in the observations, mainly as a result of information missing from the observations in the deeper atmosphere so that the dust loading at lower levels is unconstrained with large uncertainty.

The GDS onset occurs in the southern hemisphere, so the available observations during the dust storm period capture this in the assimilation at a later time than if these observations had been in the southern hemisphere. Once MCS data becomes available in the southern hemisphere again during the “cleanup” of the storm, τ ref returns toward the observed values, but again later than observed in the THEMIS data. The larger peak in the LIDO‐only assimilation suggests sedimentation in the model is not efficient enough to remove dust transported southward from the northern boundary of the GDS into the unobserved regions during the peak of the storm, a hypothesis that should be explored in future work.

Joint CIDO/LIDO Assimilation Versus THEMIS and MCS Dust Observations

Assimilating CIDO improves the dust horizontal spatial distribution, while assimilating LIDO improves the vertical distribution. By assimilating both we capture features such as the inter‐annual variability of the global mean τ ref in the THEMIS observations between GDS and non‐GDS years (Figure 4, magenta line). The global mean τ ref is reproduced during the MY28 GDS, as in the CIDO‐only assimilation (red line), and unlike the LIDO‐only assimilation (cyan line).

During the “quiet” dust season (L  = 0°–180°), jointly assimilating CIDO and LIDO gives a reasonable agreement with THEMIS observations (Figure 4, magenta line). Overestimates in the LIDO‐only assimilation during MY29 L  = 0°–90° were reduced by assimilating CIDO as well. τ ref is slightly overestimated by the joint assimilation where τ ref < 0.3 and slightly underestimated where τ ref > 2.

In the zonal‐time mean dust opacity profile compared with MCS binned observations centered at L  = 122.5° (Figure 5b, bottom panel), the joint assimilation of CIDO and LIDO produces very similar results to the assimilation of LIDO‐only (Figure 5b, fourth panel).

Comparison Between Free‐Running Model and Joint CIDO/LIDO Reanalysis

Figure 6 shows zonal mean column dust opacity τ ref in the free‐running simulation and the joint CIDO/LIDO reanalysis over MY28–29. Within each Martian year, the seasonal variability of the dust opacity in both free‐running simulation and reanalysis exhibits at least some features that are generally consistent with spacecraft observations (M. D. Smith, 2008). Global dust opacity is higher during the second half of the year, and relatively quiet during the first half of the year. In the free‐running simulation (Figure 6a), however, the active dust period in each MY lasts longer than in observations (M. D. Smith, 2009), while the dust opacity in the reanalysis (Figure 6b) shows more realistically intermittent seasonal variability within each dusty season. The peak in dust opacity is also sharper in L in the reanalysis, and tends to shut down prior to the decline in solar forcing that occurs toward the end of northern winter.

Figure 6

Seasonal evolution of the zonal mean τ ref (column integrated dust opacity (CIDO) rescaled to 610 Pa), (a) in the free‐running simulation, (b) in the joint CIDO/LIDO reanalysis. Because dust is generally well mixed in the lower atmosphere, and hence varies strongly with surface pressure, CIDO is rescaled to the 610 Pa pressure surface to account for Mars' topography.

The interannual variability in the reanalysis (Figure 6b) is essentially the same as in observations (e.g., M. D. Smith (2008, Figure 8a), M. D. Smith (2009, Figure 6) and Montabone et al. (2015)). Global dust storms (GDS) do not happen every MY, but the reanalysis successfully reproduces the observed GDS around MY28 L  ≈ 265°–310°. The initiation and duration of the GDS in the reanalysis are also consistent with THEMIS dust retrievals (M. D. Smith, 2009, Figure 6 upper panel).

Figure 8

5‐sol global mean τ ref over MY28–29 showing the free‐running simulation (green), joint column integrated dust opacity/layer‐integrated dust opacity reanalysis (magenta), and out‐of‐sample Thermal Emission Imaging System observations (triangles). Compare Figure 4 for the in‐sample observations.

The mild dusty season in MY29 following the MY28 GDS also suggests a more realistic interannual variability in the renanalysis. The free‐running simulation displays some variability, with a slightly stronger dusty season in MY28 than in MY29 (Figure 6a), but remains some considerable way away from the observations.

Verification II.: Independent, Non‐Assimilated Observations

In this section, the reanalysis from the new dust assimilation system is validated against non‐assimilated data, including the independent upward‐looking measurements from the Mars Exploration Rovers (MER) “Spirit” and “Opportunity” (Bell III et al., 2003; Lemmon et al., 2004). In order to have a more comprehensive validation, about 10% of the THEMIS and MCS data were withheld from the assimilation, and they are also used as an out‐of‐sample validation. Those withheld data were selected from 1 in every 10 of the data (for THEMIS) and of the profiles (for MCS). It would not be surprising to see that these selected THEMIS and MCS data for validation may have correlation with the data assimilated into the model, and this, to some degree, compromises their application to validate the reanalysis. The completely independent datasets from the MER landers, however, provide a complementary way of validation that does not suffer from these correlations. Hereafter, the reanalysis/assimilation used refers to the joint assimilation of CIDO and LIDO, as described in Section 4.3.

Out‐of‐Sample THEMIS Dust Observations

The distribution of non‐assimilated (out‐of‐sample) THEMIS CIDO retrievals is shown in Figure 7a. Coverage is similar to the full THEMIS data set (Figure 1a), and while it has only 10% of the data points, its distribution in latitude and L reflects the spread in the full THEMIS data set. To compare with the out‐of‐sample THEMIS data, both the observations and model results were rescaled to the 610 Pa pressure level to account for Mars' topography, and the model results were interpolated both horizontally and in time to the out‐of‐sample data points.

Figure 7

As Figure 1, but for out‐of‐sample dust retrievals.

Figure 8 shows the comparison with the global mean τ ref using out‐of‐sample observations. The free‐running simulation tracks the observations up to L  = 240° of MY28, but fails to capture the subsequent GDS. It does predict a marginally milder dust season in MY29 compared to MY28, as observed, but does not reproduce the observed dust in either case.

The reanalysis performs significantly better, capturing the MY28 GDS as well as the precursor initiation events and subsequent decay, and the interannual variability during MY29's dusty season. The magnitudes in the reanalysis are also more consistent with observations than the free‐running model, although the maximum τ ref during the MY28 GDS is still lower than observations. Measurement uncertainties in the THEMIS data may be 20% or higher, however (M. D. Smith, 2004), so the reanalysis could still be broadly consistent with the observations at the peak of the GDS.

During the “quiet” season, the free‐running model predictions generally fit the THEMIS data well, especially during MY29, at least in a global average sense. During this period both free‐running model and observations fall within the minimum observational uncertainty, which is 0.104 for the visible extinction opacity (M. D. Smith, 2009).

Correlations between τ ref in the out‐of‐sample THEMIS observations and the free‐running simulation and reanalysis are shown in Figure 9. As with the in‐sample observations, the free‐running simulation generally underestimates dust loading, mainly in the dusty season. The reanalysis produces significantly better correlations with out‐of‐sample THEMIS observations. τ ref in the reanalysis only slightly overestimates observations where τ ref < 0.5, and slightly underestimates them where τ ref > 2.

Figure 9

Scatter plots showing individual τ ref points comparing the out of sample Thermal Emission Imaging System observations with the free‐running model and various reanalyzes over the period shown in Figure 4. Colors show the data density as the number of points per square of side τ ref = 0.05. Red lines show the linear least square fit, with m the fitting coefficient, r 2 the coefficient of determination, and err the standard error in m. Black lines show m = 1. (a) Free‐running simulation, (b) joint column integrated dust opacity/layer‐integrated dust opacity reanalysis.

Out‐of‐Sample MCS Dust Observations

Figure 7b shows the distribution of out‐of‐sample MCS dust profiles. As with THEMIS there is a similar distribution pattern to the full data set (Figure 1b). The model results were averaged over several pseudo‐height ranges (0–10 km, 10–20 km, 20–30 km, 30–40 km, and 40–80 km), assuming a 10 km scale height and a 610 Pa surface pressure. Within each pseudo‐height range, the mean difference in dust opacity between the out‐of‐sample data and the joint CIDO/LIDO reanalysis was calculated for several latitude bands and is shown in Figure 10 for the “quiet” and “dusty” seasons.

Figure 10

Mean difference in dust opacity between the out‐of‐sample Mars Climate Sounder (MCS) observations and the free‐running simulation (green) and joint column integrated dust opacity/layer‐integrated dust opacity reanalysis (magenta). In each of the two seasons, the globe is split into latitude bands: (a) 90−50°S, (b) 50°−15°S, (c) 15°S − 15°N, (d) 15°−50°N, and (e) 50°−90°N. Gray dashed lines are the average uncertainties in the MCS observations.

During the “quiet” season (Figure 10a), in southern high latitudes (90°−50°S) the free‐running model significantly underestimates the dust opacity below 30 km, while the reanalysis reproduces the observations significantly better. In the midlatitudes of both hemispheres (50°−15°S and 15°−50°N) the free‐running simulation again underestimates the observations, although this underestimate decreases at higher altitude. The reanalysis generally falls within or close to the observational uncertainty, with the largest differences at 30–40 km. In the tropics (15°S − 15°N), the free‐running simulation generally underestimates the dust opacity. The reanalysis errors are generally larger than the MCS observational uncertainties except for 0–10 km where both observational uncertainties and probable systematic inadequacies of the dust lifting parameterizations are particularly acute. The assimilation is an improvement over the free‐running simulation for 10–20 km, 20–30 km, and 40–80 km, but overestimates the dust opacity for 0–10 km and 30–40 km, with absolute differences larger than those in the free‐running simulation. In northern high latitudes (50−90°N) the uncertainties in the MCS observations are about 50% larger than elsewhere. The free‐running model falls within the observational uncertainties at all altitudes. The renanalysis also falls within observational uncertainty except below 10 km.

Figure 10b shows the same for the dusty season (L  = 180°–360°). In general, the reanalysis agrees better with the MCS observations than the free‐running simulation. The maximum error in the free‐running simulation increases from south to north, which may be due to the difficulty in predicting frontal dust storms in northern high latitudes during dusty season.

In southern high latitudes, the free‐running simulation tends to underestimate the dust opacity above 10 km, and the reanalysis tends to overestimate the dust opacity below 30 km, but fall within the observational uncertainty above 30 km. In southern midlatitudes, the reanalysis is closer to the observations than in the tropics and northern midlatitudes, and at 10–20 km, 20–30 km, and 40–80 km is within or close to observational uncertainty. In the tropics, the free‐running simulation is similar to the northern middle latitudes, except below 10 km where it slightly overestimates the dust opacity. The reanalysis tends to underestimate the dust opacity between 10 and 30 km, and overestimate the dust opacity above 30 km. In the northern midlatitudes, the reanalysis underestimates the dust opacity below 30 km, but the differences are still smaller than the free‐running simulation. Above 30 km, the reanalysis overestimates the dust opacity with differences larger than the free‐running simulation. In northern high latitudes, the reanalysis falls either within or close to the MCS observational uncertainties, with lower differences at higher altitudes.

Independent Pancam Observations From Spirit

The reanalysis and free‐running simulation were compared with Spirit and Opportunity observations by interpolating τ ref horizontally and in pressure to the rover locations at Gusev Crater and Meridiani Planum respectively. Figure 11 shows these values at the locations of the two rovers during MY28–29.

Figure 11

τ ref from the reanalysis (magenta) and free‐running model (green) compared with (a) Spirit and (b) Opportunity Pancam observations. Each point is averaged over one sol.

Figure 11a shows τ ref at the Spirit rover site. During the relatively “quiet” dust season Spirit Pancam observations are normally below τ ref = 0.3. Although the free‐running simulation agreed well with THEMIS τ ref observations globally during this season (see Figure 8), at the Spirit landing site it generally underestimates the dust loading. τ ref only reaches ∼0.1 during the “quiet” season at this location. During the dusty season the free‐running simulation suggests an increase in τ ref at Spirit's location, but its increase does not match the increase in the observations.

Conversely, the reanalysis agreed better with the Spirit Pancam data. It captured the annual and interannual variability in the data well. During the “quiet” season, the reanalysis reproduces the magnitude and variation in τ ref at the Spirit landing site. Underestimates are mainly during MY29 L  = 60°–120°. The reanalysis captures the increase of dust loadings during MY29 L  = 140°–160°, but not the peak τ ref. The assimilated τ ref in Figure 11a tracks the opacity variations observed by Spirit reasonably well even during the limb staring interval of MCS during MY29 L  = 180°–270°. This indicates that the issues with calibration and limited vertical coverage have only had a limited impact, though possible systematic errors in the vertical distribution of dust cannot be ruled out. The MY28 GDS is reflected by an increase in τ ref to 3.5 at the Spirit landing site. The reanalysis reproduces the initiation and decay of the MY28 GDS, and also the variability of dust loading during MY29. Nevertheless, the reanalysis dust loading during the first peak in MY29 (L  ≈ 160°) still does not reach the maximum observed by Spirit. It is worth noting that, although the free‐running simulation fails to produce the observed amount of dust in both dusty seasons, it does exhibit some interannual variability.

Independent Pancam Observations From Opportunity

Figure 11b shows τ ref at the Opportunity rover site compared with the free‐running model and reanalysis. The observed dust loading at Meridiani Planum has a similar evolution to that seen at Gusev Crater, but with slightly higher values in general. Both the free‐running simulation and the reanalysis underestimate the peak of the MY28 GDS, though τ ref in the reanalysis is much closer to the measurements. The reanalysis also better reproduces the observed variability of dust loading during both dusty seasons, including the limb staring period in MY29 as mentioned for Spirit above, but both the free‐running simulation and the reanalysis underestimate the dust loading during the “quiet” season.

A similar discrepancy was also noticed in earlier studies when comparing datasets from different instruments. Montabone et al. (2015) found a systematic underestimate of dust opacity over the Opportunity landing site in Meridiani Planum in both TES and THEMIS datasets starting at the spring equinox, up to a factor ∼2 during northern summer, which may have been linked to the likely presence of clouds. Since the THEMIS data assimilated in this study falls within the same period, it is not surprising to find a similar discrepancy in our reanalysis compared with the Opportunity data. Lemmon et al. (2015) also raised problems with Opportunity 880 nm data around L  = 30°–130°. The source of this discrepancy remains an open problem.

Validation of the Mars Climate Database Against the Reanalysis

Reanalyses produced by data assimilation are important for verifying models against “reality” or “ground truth”. If the reanalysis is of sufficiently high quality, then it can be used as a surrogate for the real atmosphere when validating and verifying model output (e.g., ERA‐40, Uppala et al., 2005, for the Earth's atmosphere). We used the reanalysis described in previous sections to verify the Mars Climate Database v5.2 (MCD, Lewis et al., 1999; Millour et al., 2015) against our “real” atmosphere.

The quantities of interest produced in both the reanalysis and the MCD are surface pressure, surface temperature, air temperature, density, and zonal and meridional velocities. We did not compare dust diagnostics, because the reanalysis dust distribution is likely to have more high frequency variability than the MCD, due to the 3‐min period between which observations are assimilated. This will tend to stimulate short timescale variability, in a qualitatively similar way to how stochastic sub‐gridscale parameterizations also enhance variability in numerical models over a wide range of timescales (e.g., see Berner et al., 2017). These short time periods are not represented in the MCD, which is forced by dust fields changing over the timescale of at least a day. There are also significant differences between the way the dust is treated in the MCD and in the GCM used to construct the reanalysis, such as the GCM used here assumes a fixed dust particle radius of 1.5 μm, while the MCD uses a two‐moment scheme which retains information about the full dust particle size distribution (Madeleine et al., 2011). This deficiency of our model will be improved in the future.

For each 30° L period in the reanalysis (MY28 L  = 120° to the end of MY29) we computed monthly means and day‐to‐day variability in the same way as the MCD. First we interpolated horizontally from the MCD grid (5.625° × 3.75°) to the reanalysis grid (5° × 5°). Then we interpolated atmospheric quantities linearly in log  p to 30 fixed pressure levels spaced by 2.5 km up to 40 km pseudo‐altitude above a reference pressure of 610 Pa (assuming a scale height of 10 km), and spaced by 5 km above that.

The monthly mean for variable X at (longitude, latitude, pressure) position (i, j, k) is where t = 1…N includes all times within a 30° L period. Since the orbit is elliptical, N varies with season. The day‐to‐day variability is (Forget et al., 2015, Equation 5): where and are running means over 1 sol (t ± 0.5 sols) and 10 sols (t ± 5 sols) respectively. The day‐by‐day variability removes the diurnal cycle, along with any long‐term trend, leaving the variability associated with the day‐to‐day “weather”. For atmospheric quantities we calculated zonal‐monthly means at (latitude, pressure) points (j, k): where I is the number of longitude points above the surface. To compute zonal‐monthly means for the day‐to‐day variability, we computed the root‐mean‐square of the day‐to‐day variability along each latitude circle. Note this gives us a measure of the day‐to‐day variability at points along each latitude circle, rather than the spread of values along the latitude circle.

The MCD contains day‐to‐day variability as a function of position for each month and dust scenario, and monthly means as a function of local time of day at zero longitude. To obtain equivalent monthly means for comparison with the reanalysis, we averaged over all local times of day after interpolating each column to the required pressure levels.

Figures 12 and 13 summarize the differences between our reanalysis and the MCD during the second half of the Martian year. Because this paper focuses on dust, we concentrate on the dusty season from L  = 180–360°. The patterns of mean and variability differences between the reanalysis and MCD were generally similar in all seasons (with one exception, discussed below) and for all quantities when comparing MY28 with MY29, so these summary figures show only MY29. Full sets of figures for all comparisons between the reanalysis and MCD for all quantities over all months analyzed are included as Section S3 and Figures S6–S17 in Supporting Information S1.

Figure 12

Differences between reanalysis means and Mars Climate Database (MCD) means for the period L  = 180–360° in MY29, separated into 30 L segments. Positive means the reanalysis value is larger than the MCD value. Surface quantities are monthly means and atmospheric quantities are zonal‐monthly means. Gray shows points with missing data (where all points along a latitude circle are below the surface). Black lines show orography at intervals of 4 km between −4 km and +20 km above the geoid (negative values are dashed).

Figure 13

As Figure 12, but for day‐to‐day variability. Where the ratio is greater than 1, the reanalysis variability is larger than the Mars Climate Database variability over that period.

The surface temperature reanalysis was generally cooler than the MCD at most places and times, particularly between 45–60° latitude in each hemisphere. Exceptions were the polar regions and Tharsis, Arabia, and Elysium, which were persistently warmer in the reanalysis. The largest differences were near the edge of the polar icecap. The latitude at the edge of the polar icecap has a large day‐to‐day variability during a given month, as the edge of the polar cap moves over time, and so the day‐to‐day variability amplitude will transition from small (with ice) to large (without ice). This is present in both the reanalysis and MCD day‐to‐day variability (Figure S6 in Supporting Information S1). This also means the monthly mean difference between the MCD and the reanalysis is sensitive to the position of the edge of the polar cap, and so we see large differences between the reanalysis and the MCD in that region. There is a strong warm bias at southern polar latitudes between L  = 240°–330°, due to the permanent CO2 polar ice cap, which is present in the MCD but is not simulated in the version of the model used in the reanalysis.

Differences between the MCD and reanalysis surface pressure vary considerably with season. The surface pressure is persistently higher in the reanalysis in the summer hemisphere, for example, and lower in the reanalysis in the winter hemisphere. There are also rings in the difference maps around regions with the most extreme elevation changes on the planet, particularly Olympus Mons, Elysium Mons, and Hellas (Figure 12). The most likely reason is that the UK version of the Mars GCM (reanalysis) uses a spectral dynamical core, while the LMD Mars GCM (MCD) uses a grid point dynamical core. Quantities sensitive to surface elevation, such as surface pressure, will have large differences purely as a result of the topography being represented differently in the two models. The rings themselves are characteristic of the Gibbs phenomena that occur when a step function is spectrally decomposed, and are therefore likely to be spurious. Larger scale persistent differences between surface pressures in the MCD and the reanalysis may also have arisen because surface pressure in the former was specifically tuned to reproduce the Viking Lander measurements, in contrast to the model used for the reanalysis.

In the equatorial region the atmospheric temperature reanalysis is typically cooler than the MCD close to the surface, warmer around 100 Pa, cooler between 1 and 10 Pa, and warmer above 1 Pa. This pattern is repeated in most months. Both poles are typically warmer in the reanalysis than in the MCD, at least in the lower atmosphere, with a warm “tongue” in the difference maps extending into the stratospheric polar region in the winter hemisphere. There is significantly more day‐to‐day variability in the reanalysis near the surface at the winter pole.

Density differences are generally small (up to 0.2 in  log 10 ρ), but with some patterns. At low latitudes the density is generally lower in the reanalysis below 10 Pa, and higher in the reanalysis above this level. In the polar regions the density is nearly always lower in the reanalysis in the winter hemisphere. The day‐to‐day variability is lower in the reanalysis where the “warm tongue” appears in the air temperature maps. The horizontal striations in Figure 12 are likely to be artifacts: there are different vertical grids in the reanalysis and MCD, which are interpolated to a common pressure grid for comparison, and density covers several orders of magnitude, so differences in interpolation will be magnified.

Zonal velocities in the reanalysis are, in general, more westward than in the MCD. This means that the eastward mid‐latitude jets, the most prominent features of the monthly means (Figure S12 in Supporting Information S1), are weaker in the reanalysis. These differences can be quite large–up to 50 m s−1 in magnitude. At low altitudes near the equator, the zonal flow is more eastward in the reanalysis than in the MCD. The day‐to‐day variability is generally lower in the reanalysis in the winter hemisphere, and larger in the reanalysis in the summer hemisphere (Figure 13), with the exception of the winter pole near the surface, which has a high day‐to‐day variability in the reanalysis.

In general the meridional velocity reanalysis has a stronger upper‐level equatorial and midlatitude meridional circulation than does the MCD: flow away from the sub‐solar point is strengthened in the reanalysis between 0.1 and 1 Pa during the dusty season, compared with the MCD. Like the zonal velocity, typically the day‐to‐day variability is generally lower in the reanalysis in the winter hemisphere, and larger in the reanalysis in the summer hemisphere.

The only exception to the similar results for MY28 and MY29 was L  = 270–300°, which contains the build up to and peak of the MY28 global dust storm. Figure 14 shows the differences between the reanalysis and MCD during this period for both years. The main difference between the two years is that during MY28 the day‐to‐day variability is significantly larger in the reanalysis than in the MCD, when compared with the corresponding period during MY29. This applies to the surface temperature (particularly near the equator), and atmospheric temperatures, density, and both zonal and meridional velocities. Throughout the reanalysis sequence the day‐to‐day variability in the reanalysis is typically 1–2 times that in the MCD. This is likely due to the greater importance of shorter timescales in the reanalysis than in the simulations used to generate the MCD, such as a minute timescale in the temperature field, as that is the interval between successive calls to the data assimilation procedure, and this may be expected to stimulate variability on timescales shorter than a day. However, during the MY28 global dust storm period this difference is amplified. During MY28 there is a clear warm anomaly (and corresponding low‐density anomaly) between 1 and 100 Pa in the reanalysis compared with the MCD, compared with the corresponding period during MY29. A second major differences is that the monthly mean meridional velocity (i.e., the cross‐equatorial flow) between 0.1–1.0 Pa is stronger by almost 10–20 m s−1 in the reanalysis compared with the MCD during the MY28 global dust storm, while during the corresponding period in MY29 the difference is ±5 m s−1.

Figure 14

Reanalysis versus Mars Climate Database for the MY28 global dust storm period L  = 270–300° and for the corresponding period in MY29. Both the difference in monthly means and the ratio between day‐to‐day variability are shown. The units only apply to the mean differences.

Conclusions

The data assimilation system integral to the UK‐LMD Mars GCM described by Lewis et al. (2007) has been updated. That work assimilated temperature and CIDO, prescribing a vertical dust distribution using an empirical function of height (Conrath, 1975). The new scheme adds active dust lifting, transport, and deposition schemes, assuming a single dust particle size. It also assimilates vertically resolved dust profiles via LIDO, either instead of CIDO or in addition to it. This update has been prompted by the acquisition of vertically resolved dust profiles by MCS on board MRO (McCleese et al., 2010). The use of active dust transport by the winds in the assimilation whilst also assimilating dust measurements results in significant improvements in how the evolution of individual dust clouds is represented. An example contrasting assimilation with simpler binning and interpolation is illustrated in the Section S2 and Figures S4–S5 of Supporting Information S1.

When CIDO is assimilated by itself, the assimilation can reproduce the observed interannual variability of the dust horizontal spatial distribution, including the generation and dissipation of the MY28 GDS (Figure 4, red line). This inter‐annual variability cannot be reproduced by our free‐running model, although some degree of tuning can be done to reproduce the “quiet” (L  = 0°–180°) and “dusty” (L  = 180°–360°) periods during a single year. Even without assimilating any vertical information, the CIDO assimilation reduces systematic errors in the model's estimate of the dust vertical distribution. However, it misses detached dust layers that form during northern spring and summer (Figure 5b), which have been a challenge for Mars GCMs to reproduce.

Conversely, when LIDO is assimilated by itself the model can reproduce some features of the detached dust layers (Figure 5b, fourth panel), with reasonable interannual variability, which the free‐standing model is unable to reproduce. Rafkin (2012) discussed the difficulty of producing this detached dust layer in model simulations, especially in a relatively coarse resolution GCM. Similar detached dust layers were reproduced in mesoscale model simulations of “rocket dust storms” by Spiga et al. (2013) and have since been parameterized in a coarse resolution GCM (C. Wang et al., 2018) with some success. Nevertheless, being able to reproduce them in a reanalysis provides a valuable alternative means of investigating their observed characteristics and impact on Mars' atmospheric circulation. However, the limb‐viewing MCS does not continuously observe the lowest part of the atmosphere, where the dust concentration is generally highest except where there are detached layers. Consequently assimilating LIDO only does not reproduce the observed global average τ ref as well as when CIDO is assimilated by itself (see Figure 4).

Once combined together, the joint assimilation of CIDO and LIDO benefits from information about the total column dust opacity from the THEMIS data set, and the vertical distribution from the MCS data set. The evolution of the MY28 GDS is tracked well, and some features of the detached dust layers are also reproduced well.

The joint assimilation of CIDO and LIDO is a powerful tool that helps us to reconstruct the Martian climate as well as individual dust events. For example, it is difficult to retrieve a complete vertical dust distribution from MCS measurements during the MY28 GDS (Figure 1b). Using the joint assimilation, CIDO provides information to constrain the model where vertical profiles are sparse, so assimilation can map the four‐dimensional dust distribution during the MY28 GDS in the absence of complete observations. In future, direct CIDO measurements from downward‐looking instruments such as THEMIS may be supplemented by vertically integrated dust opacities obtained from MCS dust profiles, carefully extrapolated to the surface. Such an approach was used by Montabone et al. (2015) to derive daily dust scenarios for the past several Mars Years and subsequently assimilated for the OPENMars reanalysis (Holmes et al., 2020).

The reanalysis was validated against out‐of‐sample THEMIS and MCS dust observations, as well as upward‐looking MER Pancam observations at 880 nm. The reanalysis successfully reproduced the observed interannual and intraseasonal variability in the original THEMIS data, and generally improved the representation of the dust vertical distribution compared to free‐running simulations. In general, the free‐running model tends to underestimate τ ref, particularly during the dusty season, which can be signficant during major dust storm events. Hence assimilating dust observations serves to greatly reduce the model uncertainty during such events. Although free‐running simulations were able to simulate the pattern of dust loading during the quiet season, they failed to simulate the vertical dust distribution. This is consistent with the general observation that dust accumulates close to the ground, below the base of typical MCS dust profiles.

At the Spirit rover location the reanalysis captured the dust variability and intensity during both quiet and dusty seasons (Figure 11a). At the Opportunity rover location, the reanalysis captured the dust variability (Figure 11b) and improved the pattern of dust loading over the free‐running simulation during the dusty season. However, there were persistent underestimates of τ ref by the reanalysis, particularly during northern spring and summer. This is likely to be due to a systematic disagreement between THEMIS and Opportunity data (Montabone et al., 2015).

A systematic comparison between the reanalysis and the MCD during MY28 and MY29 exhibited a number of trends, although agreement between the reanalysis and the MCD was generally good. Day‐to‐day variabilities were typically larger in the reanalysis by a small factor (1–2 times), likely associated with the stimulation of fluctuations on timescales comparable with the interval between successive calls to the data assimilation procedure (3 min). In a similar manner to the impact of stochastic parameterizations of sub‐gridscale turbulence (e.g., Berner et al., 2017), this may lead to enhanced variability on a range of timescales, including those of synoptic (day‐to‐day) variability. We also found some uncertainty in the position of the edge of the polar caps in the MCD, since the main differences in surface temperature arise along those latitudes. The reanalysis also showed significantly warmer poles than in the MCD at most times of year, along with weaker mid‐latitude jets, and a stronger cross‐equatorial flow during the MY28 global dust storm.

The combined CIDO‐LIDO dust reanalysis as presented here significantly improves the estimation of Martian horizontal and vertical dust distributions over a free‐running model and over CIDO or LIDO assimilation alone. The reanalysis provides a solution generally consistent with the available Martian dust observations. The issue highlighted by Navarro et al. (2017) over the diurnal variation of dust vertical distributions has not been examined in detail here, though it does not seem to have produced the problems for temperature assimilation experienced by Navarro et al. (2017). This may be because our assimilation scheme explicitly “nudges” the model solution toward observations, which can effectively compensate at least partly for model biases and other limitations. This is in contrast to the EnKF method used by Navarro et al. (2017) which is more sensitive to model biases and imperfections. Despite this issue, which still needs further investigation in future work, our experience here suggests that assimilation has considerable potential as a tool for studying individual dust lifting events and for mapping Mars' three‐dimensional dust distribution over time. Elsewhere, we will report on further case studies using the scheme, including a southward‐moving regional dust storm during MY29, and the global dust storm during MY28.

Supporting Information S1

Click here for additional data file.