MARSWRF Prediction of Entry Descent Landing Profiles: Applications to Mars Exploration

In this paper we use the Mars implementation of the Planet Weather Research and Forecasting model, MarsWRF, to simulate the Entry, Descent and Landing (EDL) vertical profiles from six past missions: Pathfinder, Mars Exploration Rovers Opportunity and Spirit, Phoenix, Mars Science Laboratory Curiosity rover, and ExoMars 2016 (Schiaparelli), and compare the results with observed data. In order to investigate the sensitivity of the model predictions to the atmospheric dust distribution, MarsWRF is run with two prescribed dust scenarios. It is concluded that the MarsWRF EDL predictions can be used for guidance into the design and planning stage of future missions to the planet, as it generally captures the observed EDL profiles, although it has a tendency to underestimate the temperature and overestimate the density for heights above 15 km. This could be attributed to an incorrect representation of the observed dust loading. We have used the model to predict the EDL conditions that may be encountered by two future missions: ExoMars 2020 and Mars 2020. When run for Oxia Planum and Jezero Crater for the expected landing time, MarsWRF predicts a large sensitivity to the dust loading in particular for the horizontal wind speed above 10‐15 km with maximum differences of up to ±30 m/s for the former and ±15 m/s for the latter site. For both sites, the best time for EDL, that is, when the wind speed is generally the weakest with smaller shifts in direction, is predicted to be in the late morning and early afternoon.


Introduction
In order to deliver a lander or rover to the surface of a planet with the aim of better understanding the in situ atmospheric dynamics and surface geology, every mission has to go through the critical Entry, Descent and Landing (EDL) stage that directly determines the success of the mission (Braun & Manning, 2007;Raiszadeh & Queen, 2004;Steinfeldt et al., 2010;Li & Jiang, 2014). It is important to note that more than two thirds of the missions to Mars have ended in failure (e.g., Lele, 2014;Shuang & Yuming, 2011). Even though in most cases this was due to electronic/software issues, schedule constraints or human error and hence not directly related to uncertainties in the atmospheric conditions (e.g., Bitten et al., 2006;Sauser et al., 2009), a better understanding of the atmospheric variability/uncertainty may help to exclude some landing sites/seasons as well as open up more EDL options. The most recent example is the Schiaparelli Entry Demonstrator Module (EDM) that crashed landed in October 2016 (Aboudan et al., 2018). The EDL stage takes less than 10 min and occurs after a long journey to the planet (for Mars the trip takes at least 6 months). During the EDL, many complex processes, including the parachute deployment, airbag inflation, and ignition of thrusters, have to work perfectly and autonomously in a tightly synchronized manner to ensure a successful landing with no possible ground intervention (e.g., Li & Jiang, 2014;Pham et al., 2004;Way et al., 2013). As discussed in Braun and Manning (2007), and for an EDL on Mars, there are several engineering challenges that have to be addressed. For example, the fact that the Martian atmosphere is relatively thin (the density is roughly 1/100th of that of Earth's) means that entry vehicles tend to decelerate at much lower altitudes, which gives much less time for the subsequent EDL events, in particular for high elevation sites (which have lower air densities). Another complicating factor is the abundance of suspended dust. Dust has a strong impact on the thermal and dynamical state of the Martian atmosphere (e.g., Natarajan et al., 2015;Newman et al., 2002aNewman et al., , 2002bZurek et al., 1992, here after N15). A significant atmospheric dust loading, such as that observed during major dust storms, will cause the atmospheric circulation to inflate and lead to a significant increase in the atmospheric density at upper levels (e.g., Bougher et al., 2001). The Hadley Cell expands and strengthens, with major differences in the atmospheric structure including at the heights relevant to the critical Descent and Landing stages of the EDL (e.g., Martinez et al., 2017;Toigo et al., 2018). Together with atmospheric density variations, wind variations, in particular at lower elevations when the spacecraft speed is lower and hence it is more sensitive to changes in the wind, are also important. These factors together with navigation uncertainties and the specifications of the spacecraft explain why the actual landing location can be far off the center of the predicted landing ellipse (e.g., Bonfiglio et al., 2011;Knocke et al., 2004). At the surface, and even within a landing ellipse, highly variable rock abundance and rugged and variable terrain can be a concern (e.g., Bhardwaj et al., 2019). All these problems prevent the design of a standard EDL system for Mars missions and make it imperative to have a reliable estimate of the atmospheric profiles of density and winds to be encountered by the spacecraft for a better planning of the risky EDL stage and a reduction of the size of the target landing ellipse.
On Earth there are extensive studies of vertical wind profiles such as those performed by sounding balloons, vertical towers, aircrafts, and satellites (e.g., Koren & Kaufman, 2004;Schwartz & Benjamin, 1995;Tian & Lü, 2017), in addition to model simulations that are evaluated against observational data (e.g., Zhong & Fast, 2003). Model forecasts can also be combined with observational data to produce a reanalysis dataset, our best guess of the actual state of the atmosphere (e.g., Dee et al., 2011). While on Earth there is a good understanding of the near-surface winds and density profiles due to the large number of available measurements, on Mars there have only been to date six surface weather stations-Viking Landers 1 and 2 (Chamberlain et al., 1976), Pathfinder , Phoenix (Taylor et al., 2008), Mars Science Laboratory (MSL) Curiosity rover (Gómez-Elvira et al., 2012) and the Interior Exploration using Seismic Investigations, Geodesy and Heat Transport (InSight) lander (Banfield et al., 2019), and eight EDL profiles at different sites and seasons-Viking Landers 1 and 2 (Seiff & Kirk, 1997), Pathfinder (Magalhães, 1998), Phoenix (Dickinson, 2008), Mars Exploration Rovers (MER) Opportunity and Spirit (Withers & Murphy, 2009), MSL Curiosity rover (Holstein-Rathlou et al., 2016), and Schiaparelli EDM (Aboudan et al., 2018). This makes it even more challenging to have an accurate prediction of the atmospheric conditions, which in turn are needed to design the EDL phase as well as to determine the size of the landing ellipse. In addition, an accurate prediction of the atmospheric vertical profiles has implications for the early phases of assessment of the safety of one landing site versus another, which are used to decide where these extremely challenging missions will land and operate. With this in mind, it is useful to consider the key EDL phases and how they relate to atmospheric conditions. The EDL comprises the EDL steps. The Entry phase begins when the vehicle reaches the top of the Martian atmosphere (at roughly 120-130 km above ground level, AGL) and ends with the deployment of the supersonic parachute, lasting about 4 min. During this phase the speed of the spacecraft is reduced from about 4,000-7,000 m/s to about 400 m/s through friction with the atmosphere. The parachute deployment takes place once a specified threshold (such as Mach number and altitude) is met, typically around 10 km AGL. After the vehicle quickly decelerates from supersonic (~400 m/s) to subsonic (~80 m/s) conditions, the heatshield is jettisoned, at about 6-7 km AGL, and an accurate measurement of the vehicle's altitude and velocity relative to the surface is taken before proceeding with the powered descent. The timing of the jettison of the heatshield depends on two constraints: clear separation from the flight system to prevent recontact and sufficient distance (at least 15 m) so as not to block more than one of the beams of the terminal descent sensor. The powered descent phase begins with the backshell separation at about 2-3 km AGL and comprises a series of events separated by altitude-velocity waypoints to ensure a controlled operation. Prior to hitting the surface, the thrusters are ignited to slow the vehicle's speed of descent and the airbags are inflated to cushion the surface impact. Only when it comes to rest the airbags are deflated and retracted and the lander or rover starts operating. It is important to stress that the architecture of the different EDL stages has evolved from the Viking to the InSight missions (e.g., San Martin et al., 2013). For example, and regarding the landing stage, the Phoenix, MSL, and InSight missions, due to their weight, have not used an airbag-assisted landing. As landers, Phoenix and InSight used thrusters attached to them (e.g., Abilleira et al., 2014;Prince et al., 2011), while the MSL, in order not to be burdened by thrusters on the surface, used the "sky crane" touchdown system (e.g., Steltzner et al., 2006). What is more, new technologies are likely to be used in future missions to Mars (e.g., Jiang et al., 2013). A detailed summary of the processes involved in the EDL is given in Pham et al. (2004), Braun and Manning (2007), Prince et al. (2011), Way et al. (2013, and Li and Jiang (2014).
The initial EDL operations take place at very high spacecraft speed, and hence, the wind and density variations have a limited impact. However, they are particularly important at the moment of deployment and operation of the final parachutes, as well as when the retro-propulsion rockets are activated or when parts of the landing structure are jettisoned to be driven far away from the main lander. During these final stages of descent and landing, the speed of the spacecraft is low and the uncertainties in the ratio of atmospheric wind speed to the spacecraft speed become more relevant. Strong wind speeds and/or significant shifts in the wind direction may interfere with these events, which normally take place in the lowest 10 km above the surface. For instance, for the very last stage of the MSL landing phase, the constant deceleration segment, the Powered Descent Vehicle was decelerated from 20 to 0.75 m/s in a controlled way (Prakash et al., 2008). In addition to the strength of the wind, fluctuations in the wind direction can also affect the spacecraft attitude stability control system or the safety of certain maneuvers like heatshield jettisoning. The parachute is generally deployed at a distance of about 10 km from the surface (e.g., N15; Holstein-Rathlou et al., 2016), and thus, special attention must be paid to the lower layers of the atmosphere. As a result, low horizontal wind speeds with small shifts in wind direction are preferred during the EDL stage and will be considered in this work as the most favorable atmospheric conditions for EDL.
Despite all the challenges inherent to the EDL, it provides a once-per-mission chance to retrieve vertical profiles of variables such as density, pressure, and temperature at very high vertical and horizontal resolution, especially in the lower atmosphere, making it scientifically valuable even in comparison with those estimated from orbital assets such as the Mars Global Surveyor (MGS) Thermal Emission Spectrometer (TES) and the Mars Reconnaissance Orbiter Mars Climate Sounder (Holstein-Rathlou et al., 2016). As explained in detail in Aboudan et al. (2018), the atmospheric acceleration that is measured during the descent is directly related to the atmospheric density through the drag equation, which allows its retrieval. Given the density profile, the pressure can be derived from the hydrostatic balance. The temperature profile is then computed using the ideal gas law: p = ρRT where p is the pressure (Pa), ρ is the density (kg/m 3 ), R is the specific gas constant (for the Martian atmosphere it is about 192 J kg -1 K -1 ), and T is the temperature (K). In addition, the horizontal wind speed and direction can also be derived from the accelerometer data recorded during the EDL stage (Aboudan et al., 2018) or by examining the impact points of the jettisoned hardware (Paton et al., 2018), although some assumptions have to be made.
In addition to providing insight into the vertical structure of the atmosphere, the reconstructed atmospheric profiles can be used for evaluation of the performance of numerical models. This has been done for the Phoenix Mars Lander (e.g., Withers & Catling, 2010), MER Spirit and Opportunity (e.g., N15), MSL (e.g., Holstein-Rathlou et al., 2016), and Schiaparelli EDM (Aboudan et al., 2018) missions. Numerical simulations have also been conducted in support of planned Martian missions (e.g., Kauhanen et al., 2008;Michaels & Rafkin, 2008;Rafkin & Michaels, 2003;Rafkin et al., 2004;Tyler et al., 2008). Setting up and evaluating atmospheric simulation models is important for the engineering design of EDL phases because once the atmospheric plausible range of realistic, and extreme, scenarios are given, engineers can implement Monte-Carlo modeling of different entry paths through random or selected mesoscale profiles and add high-frequency turbulence to simulate how the automatic control system will perform. Finally, it is also important to predict the landing site ellipse and to calculate the probabilities of landing at locations where mobility would be restricted, for example, by topographic obstacles and rock abundance, or where horizontal landing may be compromised.
In this paper the MarsWRF model is run for the landing sites of six past (namely, Mars Pathfinder, MER Spirit and Opportunity, Phoenix Mars Lander, MSL, and ExoMars 2016 Schiaparelli) and two future (ExoMars 2020 surface landed platform and rover and Mars rover 2020) missions to the planet. The InSight mission is not considered as the EDL data have not been made publicly available at the time of this study. The main goal is to investigate the quality of its predictions of the EDL profiles and their sensitivity to the vertical distribution of dust. The focus on dust is justified given its very important role in the vertical structure of the atmosphere as shown, for example, by N15, also using the MarsWRF model. The implications of the model predictions for the EDL of upcoming missions to Mars will also be discussed.
New ambitious plans are being set up by the National Aeronautics and Space Administration (NASA) and the European Space Agency (ESA) for the near future, first regarding Mars Sample Return (MSR; e.g., McLennan et al., 2012;Stefaan & Beaty, 2018) and then, eventually, the human exploration of Mars (Grant et al., 2018). For MSR a new technology, which has never been demonstrated, needs to be put in place, namely, the ascent of a landed spacecraft from a planet (i.e., an inverse EDL) where the initial winds may also affect the launch from Mars. Before sending humans to Mars there needs to be a number of spacecraft sent to the planet demonstrating a successful execution of both the EDL phase and its inverse for spacecraft of increasing mass and complexity. Furthermore, the distance between the actual landing site and the center of the landing ellipse has to progressively be reduced: for example, Spirit and Opportunity, who landed during the dust storm seasons, overshot the center of their landing ellipses by about 10.1 and 24.6 km, respectively, whereas with more favorable atmospheric conditions, MSL overshot the center of its landing ellipse by about 2.4 km (Cantor et al., 2019). This further motivates a dedicated study of the different atmospheric profiles at landing sites on Mars to facilitate the future engineering predesign phase of these missions, which in most cases may be led by space agencies (national or international) but also eventually by private companies (such as, e.g., SpaceX) or joint private/public missions. This manuscript is divided into five sections. Section 2 provides details about the model configuration used. The results of the MarsWRF simulations for past Mars missions are shown in section 3, while in section 4 the focus is on future missions to the planet, in particular ExoMars 2020 and Mars 2020. A discussion of the main findings is presented in section 5.

Experimental Setup
The Martian implementation of the Planet Weather Research and Forecasting model (Fonseca et al., 2018;Newman et al., 2017;Richardson et al., 2007), MarsWRF, is used to simulate the atmospheric conditions at the time of EDL of six past and two planned missions to Mars. The topography, surface albedo, and thermal inertia used in the model simulations are derived from measurements taken by instruments on the MGS spacecraft. The topography is interpolated from the 1/64°global topographic grid obtained from Mars Orbiter Laser Altimeter (MOLA) estimates (Smith et al., 2001). The albedo map employed is that for Martian Year 26 (MY26; Clancy's et al., 2000, convention for MYs is used here) derived from measurements by TES. For thermal inertia the TES apparent nightside thermal inertia map is used (Putzig & Mellon, 2007). A simple CO 2 microphysics scheme that accounts for the sublimation and deposition of CO 2 and makes the assumption that all condensed CO 2 is directly deposited onto the surface is employed in the model. The total CO 2 inventory and polar cap properties are adjusted so as to produce a pressure cycle that best matches the Viking Lander 1 and 2 observed pressure curves (Guo et al., 2009). Thermal conduction within the subsurface is represented by a 15-layer implicit scheme. Following Richardson and Wilson (2002), MarsWRF is not set up with a permanent CO 2 ice cap at the southern pole, but its effects are accounted for by keeping the ground temperatures at latitudes below 82.5°S equal to the CO 2 frost temperature. The planetary boundary layer (PBL) scheme used is the Medium Range Forecast (MRF; Hong & Pan, 1996) PBL scheme, named after the model in which it was first used, which includes both local and nonlocal mixing within the PBL and only local mixing above the boundary layer. The radiation scheme employed in the model runs is a k-distribution radiative transfer model described in Mischna et al. (2012) with a non-local thermodynamic equilibrium adjustment to the long-wave radiation flux based on López-Valverde et al. (2000). Further details regarding the model's grid, dynamical core, and physics parameterization schemes can be found in Newman et al. (2017) and Fonseca et al. (2018).
The two main goals of this study are (i) to investigate the sensitivity of the atmospheric profiles to changes in dust loading at different locations on the planet and for different times of the year, in particular for the six past missions for which the reconstructed profiles are available, and (ii) to use MarsWRF to predict the EDL conditions that may be encountered by the two upcoming missions to Mars. If one assumes zero knowledge of the actual atmospheric dust loading at the EDL time, the most representative dust scenarios are "average" scenarios, that is, those that are not specific to a given MY. In order to study the sensitivity to the dust loading, two yearly repeatable global dust distributions are used in this work: the MCD MGS dust scenario ("MCD"; Montmessin et al., 2004;Toigo et al., 2012) and one derived from TES limb opacities ("TES"; . The MCD dust scenario aims to fit most of the thermal profiles observed during the MGS mission representing our current "best guess" of the dust on Mars as observed by the MGS. It is constructed based on the assumption that the dust mixing ratio is constant from the surface to a (spatially and temporally varying) height, above which it rapidly declines. The TES dust scenario is obtained by applying a full trilinear interpolation (with linear extrapolation) to the TES limb opacities averaged over L s 110°MY24 -L s 80°M Y27, together with a sinusoidal time of day fit with a peak and trough at TES day and night measurements, respectively. In order to construct this scenario, the TES limb opacities for the roughly 3 MY period are first extracted, with the total column opacity at a given longitude, latitude, and time normalized to match the concurrently observed TES nadir observations. The median dust optical depth for the 3 MY period is then estimated to generate the "annual average" dust scenario. It is important to note that the TES-climatology scenario considered here is not the same as that used in N15. In particular, in the adjustment done to the limb opacities so that the column optical depth matches the TES nadir observations, the scaling in the TES-climatology scenario used here is applied at all vertical levels, whereas in N15 it is only performed near the surface where the uncertainties in the opacity measurements are higher. Figures 1a and 1c show the seasonal and latitudinal distribution of the dust optical depth normalized to 7 hPa, used as proxy for atmospheric opacity in the model, for the MCD and TES dust scenarios. In both the dust opacity increases during the dust storm season from L s 180°to 270°and decreases afterward with a clearer atmosphere around the boreal summer solstice (L s 90°). Figures 1b and 1d show the height at which the dust optical depth decreases by a factor of e with respect to its value in the lowest model layer for the two dust scenarios. This height is estimated in the following way: (i) for the global grid of each simulation, and at a given sol (i.e., L s value), latitude, and longitude, the height at which the dust optical depth drops by a factor of e from its value at the bottom layer is estimated, and (ii) for a given latitude and L s , the average height over all longitudes is plotted. As can be seen, the dust optical depth generally decreases faster with height in the MCD dust scenario, except mainly in the high-latitude regions in the local summer season. In other words, at higher elevations the TES dust scenario typically features increased dust loadings, as also noted by N15. Hence, these two scenarios have not only different column opacities but also vertical distributions of dust, as they are designed to match different sets of observations: temperature profiles, for the MCD dust scenario, and the TES nadir column opacities, for the TES dust scenario. As stressed in N15, both the total column and vertical distribution of dust are crucial in determining the atmospheric thermal structure. In Figure 2 the dust optical depth profiles for the two dust scenarios around the EDL time of each of the six past Mars missions considered in this work are given and will be used for discussion in section 3. They are consistent with the plots given in Figure 1 and with the description presented above.
While a complete study would have required additional dust scenarios, not just those tailored for particular MYs like the TES MY26 used by N15 but also extreme ones, for example, representative of regional or even global dust storms to investigate their impact on the EDL profiles, the choice of these two "average" and "storm-free" scenarios is justifiable as (i) it allows for an investigation of the sensitivity of the atmospheric response to changes in the dust distribution at different locations on the planet and times of the year, and (ii) the EDL of ExoMars 2020 and Mars 2020 will take place outside the dust storm season (L s 19°).  For each experiment, the model is set up in a five-nest configuration centered on each mission's landing site shown in Figure 3. The outermost grid is global at a spatial resolution of 2°× 2°(roughly 118.33 km) with the other four domains having resolutions of about 39.44, 13.15, 4.38, and 1.46 km, respectively. Due to the decreasing zonal grid spacing with latitude, and in order to avoid the use of very small time-steps, in the outermost domain and for latitudes polewards of 70°, a zonal filter that damps high frequency waves is employed (Richardson et al., 2007). Given this, for the Phoenix EDL simulations this grid is rotated so that the North Pole sits on the computational equator. In all experiments the global grid is first run for two MYs with the first MY regarded as spin-up and the output discarded. Restart files from the second MY are kept every 10 Martian solar days (sols) and subsequently used to initialize the nested runs. As in Fonseca et al. (2018), the global grid is run hydrostatically, while the other four are run in nonhydrostatic mode. The time steps used in the five grids are 1 min, 20 s, 10/3 s, 10/9 s, and 10/27 s, respectively, with the output stored every 1 hr for the first three, every 30 min for the fourth grid and every 5 min for the fifth grid. This last grid is run for eight sols with the output of the first discarded as spin-up and that of the other seven used for analysis. Regarding the computational resources, the global simulations take about 7 days to complete with 48 processors, while the nested runs finish within 4 days with 144 processors on the High Performance Computing Center North (HPC2N) Abisko cluster. In the vertical, 45 levels are considered with the model top at 120 km. The approximate height AGL of the center of the vertical layers is given in Table 1. In all grids, and in the top 4 km, Rayleigh damping is applied to the wind components and potential temperature on a time-scale of 5 s (Skamarock et al., 2008). In particular, the zonal mean of the horizontal wind components and potential temperature are damped in the three highest model levels on time-scales of 2, 6, and 18 days, respectively, to zero for the wind components and to 140 K for the temperature. Given this, the model output at the top level is discarded, and hence, the model predictions are only shown from the surface to~35 km AGL. It is important to note that this set of vertical levels can be used for the purpose of this work. As explained in the previous section, the atmospheric conditions in the lowest 30 km are more important for EDL operations given the slower speed of the spacecraft (e.g., Way et al., 2013).
The EDL profiles of temperature, density, and pressure for the Pathfinder, MER (Spirit and Opportunity), Phoenix, and MSL missions are available at the NASA's Planetary Atmospheres Node of the Planetary Data System website (PDS; http://pds-atmospheres.nmsu.edu/). The EDL profiles for the Viking Lander 1 and 2 missions are also available on the PDS website, but these data are not certified and so are not considered here. The reconstructed EDL profiles for the Schiaparelli EDM are shown in Aboudan et al. (2018) and provided by the first and corresponding author of that paper. As pointed out by Withers and Smith (2006) and N15, care should be taken when using the estimated EDL profiles for model evaluation. Factors such as the estimation of the atmospheric wind speed, angle of attack, and drag coefficient of the spacecraft are known to cause errors in the reconstructed profiles, in particular for elevations below 15 km. As a result, whenever available, the uncertainty estimates will be considered, and both a qualitative and quantitative assessment of the model predictions will be performed.

EDL Simulations for Past Mars Missions
In this section a discussion of the results of the MarsWRF experiments for six past missions to the Martian surface is presented. For each landing site the model is run with two prescribed "average" dust scenarios, MCD and TES, to investigate the sensitivity of the vertical profiles of density, pressure, and temperature (and horizontal wind components for the Schiaparelli EDL for which the reconstructed profiles are available) to the distribution of dust. The MCD and TES dust optical depth vertical profiles for each simulation are given in Figure 2.

Mars Pathfinder
More than 20 years after the Viking missions, the first spacecrafts to successfully operate on the surface of Mars, Pathfinder was launched in December 1996 touching down on 4 July 1997 . The landing site was Ares Vallis (19.33°N, 33.55°W), located in Chryse Planitia, and the EDL took place in the local midsummer season (L s 143°) around 03:00 Local Mean Solar Time (LMST). Having lasted for nearly three times the planned lifetime of 30 sols, this mission helped to better understand the atmospheric conditions at Ares Vallis as well as the geomorphology and mineralogy of the site (e.g., Bell et al., 2000;Hviid et al., 1997;Morris et al., 2000;Smith et al., 1997). Figure 4a shows the observed (black) and model-predicted (red for the MCD and blue for the TES dust scenario) vertical profiles of density, pressure, and temperature from roughly 9 to 35 km above the landing site. For the observations, shown is the mean ±3σ uncertainty (Magalhães, 1998), whereas the MarsWRF curves show the range of values (i.e., minimum to maximum) obtained for the seven MY sols 299 to 305 (L s~1 43°) from 02:00 to 04:00 LMST (i.e., the 2-hr period centered on the landing time).
For both the MCD and TES simulations, the model-predicted density decreases slightly slower with height compared to observations below 12.5 km but faster above 15 km. As far as the temperature vertical profile is concerned, from about 9 to 16 km the predicted and observed temperatures are generally within 10 K of each other but have opposite trends: while the latter increases with height, indicating the presence of an inversion, the MarsWRF temperature decreases at an average lapse rate of 2.5 K/km, the typical lapse rate of the Martian atmosphere. According to Magalhães et al. (1999) and Colaprete and Toon (2000), the steep increase in temperature between roughly 10 and 16 km may be due to radiative cooling associated with the presence of water ice clouds that are not represented in the model. This is not surprising, as at this time of the year the equatorial or aphelion cloud belt is present mostly in the region between about 10°S and 30°N which includes Ares Vallis (e.g., Smith, 2004;Wang & Ingerssol, 2002). From about 16 km to the highest vertical level considered,~35 km, in both MarsWRF and observations, the temperature decreases with height. At these heights the amount of suspended dust has an important impact on the thermal structure of the atmosphere (Haberle et al., 1997; Magalhães   , 1999). The decreasing temperatures suggest a gradual reduction in the concentration of dust with the model being systematically colder by up to~20 K (the MCD experiment gives slightly better results). The lower temperatures and consequently higher densities predicted by MarsWRF indicate a clearer atmosphere in the model compared to that observed. A comparison of Figure 4a with Figure 2 suggests that the warmer profile in the MCD simulation arises from the generally higher atmospheric dust loadings compared to the TES run. It is important to note that the general good agreement between the observed and modeled pressure profiles results from compensating discrepancies in the temperature and density profiles: below~12.5 km the former is overpredicted and the latter underpredicted, while the opposite is true above~15 km.

MER: Spirit and Opportunity
Roughly 18 years after the Pathfinder mission, NASA sent two robotic rovers to opposite sides of Mars: Spirit, which landed at Gusev Crater (14.57°S, 175.48°E) on 4 January 2004 (L s 328°), and Opportunity, which landed at Eagle Crater in Meridiani Planum (1.95°S, 5.53°W) on 25 January 2004 (L s 339°). Having operated for more than 2200 sols, Spirit detected dust devils, monitored the dust opacity and the lower boundary layer at Gusev crater, and allowed for a better understanding of the geology of the site where soils are dominated by olivine-bearing basalts with evidence of an aqueous environment in the past Greeley et al., 2006). Liquid water was also likely to have been present in the past at Meridiani Planum where Opportunity landed and operated for more than 5,000 sols until June 2018, but the current conditions are very different with a predominantly arid, acidic, and oxidizing environment . Both missions are regarded as a tremendous success (Erickson et al., 2007). Figure 4b shows the model predictions for the MCD and TES dust scenarios and observed density, pressure, and temperature profiles for the Spirit EDL. The observed curve shows the mean values ±1σ uncertainty (Withers & Murphy, 2009), whereas the model curves indicate the range of values from 13:00 to 15:00 LMST (the EDL took place at 14 LMST) for the seven MY sols 606-612 (L s~3 28°). As opposed to the Pathfinder case, the observed temperature changes very little below 30 km and decreases slowly with height above that level to about 35 km. This indicates a relatively dusty atmosphere, which may be associated with the moderate dust storm that occurred on Mars shortly before the atmospheric entry of this mission (Withers & Smith, 2006;N15). As the atmosphere in MarsWRF is less dusty compared to that observed, the differences between the predicted and observed temperatures are rather large in this case, in particular above~15 km. The temperature bias exceeds 50 K in the MCD simulation, with the densities also higher than those observed. Given the higher amounts of dust loading and column opacity in the TES simulation (Figure 2), it generally gives a slightly better agreement with the observed profiles. As with the Pathfinder simulations, the MarsWRF biases in temperature and density to a certain extent offset each other, with the model giving a good simulation of the observed pressure profile. Figure 4b but for the EDL of the Opportunity mission. As for the Spirit case, in the observed data, the mean values are shown with a ±1σ uncertainty (Withers & Murphy, 2009). The touchdown was at 13:20 LMST and so the range of values shown for the model data are those obtained from 12:20 to 14:20 LMST for the seven MY sols 625-631 (L s~3 39°). In contrast with the Spirit profile, however, the temperatures decrease steadily with height from~16 to 35 km (except for a weak inversion from about 23 to 25 km). This is consistent with the observed lower local and global-scale atmospheric dust content for the Opportunity's EDL compared to the Spirit's EDL (Withers & Smith, 2006;N15). The magnitude of the model biases does not exceed about 20 K for the TES but can reach 30 K in the MCD simulation, in line with the dust opacity profiles given in Figure 2. Two notable features in the observed temperature profile below 16 km are not simulated by MarsWRF: a nearly isothermal layer from about 12 to 16 km and a strong inversion from 8 to 12 km. The presence of radiatively active water ice clouds and their potential coupling with thermal tides may explain that inversion (Hinson & Wilson, 2004), although it may also be an artifact of the retrieval technique (Withers & Smith, 2006). The nearly isothermal layer and the weaker inversion between 23 and 25 km may be associated with a higher dust abundance not represented in the model when forced with two "averaged" prescribed dust scenarios.

Figure 4c is as
N15 also investigated the sensitivity of the MarsWRF predicted temperature profile to the atmospheric dust loading for the Spirit and Opportunity EDLs. The results they obtained with the MCD dust scenario are similar to those reported here, but for the TES run the model temperature bias shown in Figures 4b and  4c is of a larger magnitude compared to that obtained by N15, roughly twice as large. A possible explanation is the different ways in which the TES-climatology dust scenario is constructed, highlighted in section 2. In addition, N15 considered the TES MY26 dust scenario, which corresponds to the MY when the EDL of these two missions took place. They found that with this dust scenario the model predictions are in closer agreement with the observations, stressing the importance of both the dust vertical distribution and column opacity for the simulation of the EDL profiles. In addition, and as highlighted, for example, in Magalhães et al. (1999), other factors, such as the presence of water ice clouds not captured by MarsWRF as it is run without a water cycle, can also account for some of the differences between the observed and predicted temperature profiles.
The model response obtained with the MCD dust scenario for the Spirit and Opportunity EDLs shown in Figures 4b and 4c are comparable to those obtained by Montabone et al. (2006) with their "baseline" dust scenario, similar to the MCD dust scenario considered here, with a clear cold bias for heights above~15 km. However, when Montabone et al. (2006) used a dustier scenario based on the Viking Lander dust opacity observations with the large dust storms removed, a much better agreement with the observed profiles is seen. This further stresses that the cold bias in MarsWRF above~15 km may indeed be due to a reduced dust loading in the model.

Mars Phoenix Lander
After the failure of the Mars Polar Lander mission in 1999, the Phoenix mission was the first to successfully land in the Martian polar region. The landing site was in Green Valley, at about 68.22°N and 125.75°W, in a region known as Vastitas Borealis. As was the case for the other missions, Phoenix exceeded its design lifetime of 90 sols by about 67 sols and achieved the proposed scientific objectives of allowing for a better understanding of the Martian North Polar environment (e.g., Arvidson et al., 2009;Cull et al., 2010;Ellehoj et al., 2010). Figure 5a shows in black the observed vertical profiles of density, pressure, and temperature from about 14 to 35 km above the landing site with a ±1σ uncertainty (Dickinson, 2008). In red and blue the MarsWRF predictions for the MCD and TES dust scenarios are given, respectively. The touchdown took place in late boreal spring (L s 76°) at around 16:30 LMST. As a result, the model curves show the range of values from 15:30 to 17:30 LMST, a 2-hr period centered on the landing time, for the seven MY sols 164-170 (L s~7 6°). Observations and numerical simulations indicate a well-mixed boundary layer up to 4 km in the afternoon and a strong surface inversion at night with potential radiation fog and water ice clouds at 4 km AGL (Savijärvi & Määttänen, 2010;Whiteway et al., 2009). These features are not seen in the EDL profile, as data are only available for heights above 14 km and would also not be predicted by MarsWRF with its present configuration. For the heights shown in Figure 5a the amount of dust present plays the largest role in the vertical structure of the atmosphere Magalhães et al., 1999. As was the case for the other missions considered, with both dust scenarios, the atmosphere is clearer in the model compared to that observed, which is reflected in the predicted colder temperatures and higher atmospheric densities. The reduced temperature biases in the MCD simulation are consistent with the larger opacities in the model, as seen in Figure 2.

Mars Science Laboratory
The MSL is the first mission to land on a site with significant local topography on Mars. It touched down on 6 August 2012 (L s 151°) at about 15 LMST in the north-western part of Gale Crater (4.49°S, 137.42°E) and delivered the rover Curiosity that has been operating ever since. Gale is a 154-km-wide crater with a large central mound estimated to be~3.6-3.8 billion years old (e.g., Anderson et al., 2018;Thompson et al., 2011). Measurements by Curiosity have helped to better understand local, regional, and infer on the global features of the Martian atmospheric dynamics (e.g., Fonseca et al., 2018;Haberle et al., 2014;Martín-Torres et al., 2015;Newman et al., 2017;Rafkin et al., 2016). Figure 5b shows the MarsWRF predicted profiles for the seven MY sols 316-322 (L s~1 51°) and 2-hourly period centered on the landing time (here 14:00 to 16:00 LMST). The reconstructed profiles of density, pressure, and temperature from about 12 to 35 km above the landing site with a ±1σ uncertainty (Holstein-Rathlou et al., 2016) are shown in black. For the range of elevations considered (~12-35 km), the model tends to underestimate the pressure and the temperature (the latter by about 10-20 K), with these errors partially compensating and leading to a better agreement between the observed and modeled density profiles. The dust opacity vertical profiles for the MCD and TES distributions for this EDL are the closest for all six past missions considered (Figure 2), which is consistent with the comparable temperature, pressure, and density profiles for the two simulations given in Figure 5b. The weak inversion in the reconstructed temperature profile around 30 km may be associated with the thermal tides and their modulation by the atmospheric dust content (e.g., Lee et al., 2009). The lower vertical resolution of the model grid at those heights and the use of prescribed "averaged" dust distributions may explain its absence in the MarsWRF profile. The "kinks" in the reconstructed profiles around 14.5 and 18 km are likely associated with EDL events in preparation of the parachute deployment, including the ejection of the Entry Ballast Masses and consequent change in the angle of attack (Holstein-Rathlou et al., 2016;Karlgaard et al., 2014).

Exomars 2016: Schiaparelli
As the first mission of the ExoMars program, the Trace Gas Orbiter and the Schiaparelli EDM were launched in March 2016. While the Trace Gas Orbiter was successfully placed into orbit around the planet, following a successful entry and descent phase the Schiaparelli EDM crashed on the surface of Mars on 16 October 2016 (L s 245°). The crash has been attributed to a failure of the Guidance and Navigational Control subsystem that led to an early release of the back-shell and the parachute and a premature deactivation of the thrusters (Tolker-Nielsen, 2017). In any case, and despite the crash landing, the data transmitted in real time by the Schiaparelli EDM during the entry and descent stages has been used to reconstruct the vertical profiles of density, pressure, temperature, and horizontal winds (Aboudan et al., 2018) used here for evaluation. Figure 6 shows the reconstructed and model-predicted vertical profiles of density, pressure, and temperature, as well as horizontal wind components, speed, and direction. The red and blue curves give the range of values predicted by the model for the seven MY sols 473 -479 (L s~2 45°) for the MCD and TES dust scenarios, respectively, from 13:00 to 15:00 LMST (the landing time was around 14 LMST). The vertical coordinate used here is not AGL but Above MOLA Radius (AMR), the altitude with respect to the constant MOLA radius of 3396 km. The AMR elevation is given by the AGL elevation adding the terrain height for the closest model grid-point to Schiaparelli's crash site, which is about -1.45 km. It is important to note that the uncertainty for the density and temperature between 4.5 and 7.5 km AMR can be rather large, at times exceeding 200 K for the latter. This indicates that the estimates in this range are not reliable. In order to filter out those measurements, for all vertical levels and fields for which the uncertainty exceeds 25% of the mean value, the reconstructed estimates are not plotted. No uncertainty is provided for the horizontal wind estimates. Figure 6a shows a good agreement between the reconstructed and model-predicted EDL profiles in particular for the TES experiment and above~18 km. Below this level, the temperature profile in the TES run is colder by as much as 30 K, while in the MCD run the biases are roughly 10 K lower. The difference in the two model responses can be understood by looking at the dust opacity profiles given in Figure 2. Above~12 km, the dust opacity in the TES simulation decreases slower with height compared to that in the MCD run, with the larger amounts of dust in the former consistent with the considerably higher atmospheric temperatures. It is interesting to note that the dust opacity values for this EDL case are higher throughout the whole column in comparison with the other five missions considered in this work. In fact, this is the only case for which MarsWRF is successful in capturing the reconstructed temperature profile above~18 km, stressing that its underprediction in the other cases may indeed be related to reduce amounts of dust in the atmosphere. On the other hand, below~12km, the dust loading is higher in the MCD simulation, in particular at lower levels, which is consistent with the warmer atmospheric temperatures. The colder temperatures in the TES run below~12 km may arise from the increased dust loading overhead, and the vertical distribution of dust in this simulation may also explain the relative minimum in the temperature profile at~11 km, seen in the observed profile at roughly the same level but with a larger amplitude. In the reconstructed profiles, this sudden change in temperature is accompanied by a change in density, not seen in the MarsWRF plots, which Aboudan et al. (2018) attribute to an atmospheric feature that is not captured by the model. The density profiles below~12 km and pressure profiles below~5 km are generally well simulated, but above those height they are both lower than those observed.
For this mission, the EDL profiles of the horizontal wind components for heights between~2.5 and 9 km are also available for comparison. The estimated profiles based on the observed data and MarsWRF-predicted wind profiles are shown in Figure 6b. Despite not being able to simulate the fine-scale structure present in the observed profiles (e.g., the relative maxima at about 3.8, 5.4, and 7.6 km), which may be an artifact of the reconstruction procedure, the predicted maximum wind speeds are the same as in the observations, about 20 m/s. This result highlights the value of the MarsWRF predictions when forced with two dust scenarios effectively averaged over multiple MYs, and hence that do not include year-specific details such as regional dust storms, for a variable that is very important to the EDL stage of a mission. The wind tends to blow from an easterly direction, except mostly between 4 and 4.5 km, below 3 km and around 8.5 km AMR. MarsWRF generally captures the prevailing wind direction well, even though below 4 km it shows a large range of wind directions.

EDL Simulations for Upcoming Mars Missions
In the previous section it is concluded that for six of the past missions to Mars, MarsWRF has some skill in simulating the observed EDL profiles. The main bias is a tendency for colder temperatures and higher shows the vertical profiles of the zonal, meridional and total horizontal wind speed (m/s) and horizontal wind direction (º) from 2.5 to 9 km AMR. No uncertainty is available for the reconstructed horizontal wind components. In all panels, and for the model predictions, the range of values obtained in a 7-sol simulation around the EDL time and in a 2-hourly period centred on the landing time is plotted. In the pressure and density plots, the horizontal axis is in logarithmic scale. densities in particular for heights above 15 km, possibly due to a clearer atmosphere. In this section the model predictions for the vertical profiles of density and horizontal wind, two crucial variables that have to be well simulated for the success of the EDL of a mission (e.g., Knocke et al., 2004), for the candidate landing sites of the ExoMars 2020 (Oxia Planum) and Mars 2020 (Jezero Crater) missions are discussed.

ExoMars 2020 (Oxia Planum)
The proposed landing site for ExoMars 2020 is Oxia Planum (18.20°N, 335.45°E), a region that is particularly interesting from a geological point of view comprising both younger volcanic and older sedimentary rocks (Pajola et al., 2017). The landing is expected to occur in the boreal spring season (L s 19°).  above the surface where the wind speed is at a minimum. In order to better highlight the differences between the two experiments, in Figure 7b the ratio of the densities and the horizontal wind speed differences are given. Here the range of values obtained in the seven sols is plotted. While in the lowest 15 km the densities produced by the two simulations are comparable, above~20 km they are smaller in the MCD experiment with maximum differences of up to 10%. The density profiles are consistent with the dust opacity profiles shown in Figure 7a, as the atmosphere is dustier in the MCD experiment. For heights of up to~8 km the wind speed differences between the two runs are generally within ±8 m/s but between 15 and 30 km, and except at 18 LMST, they can reach ±30 m/s. This indicates a strong sensitivity of the horizontal wind speed to the atmospheric dust loading and has important implications for the planning of the EDL phase of ExoMars 2020. When taken as a whole, and according to the model predictions, the EDL of ExoMars 2020 is likely to be rather challenging, with more favorable conditions, that is, when the wind speed is the weakest with smaller shifts in wind direction, in the late morning and early afternoon.

Mars 2020 (Jezero Crater)
The proposed landing site for Mars 2020 mission is Jezero Crater (18.4463°N, 77.4565°E), with the landing also expected around L s 19°. This crater has been found to be both geologically and astrobiologically compelling as it likely hosted a stable body of standing water in the past (e.g., Goudge et al., 2017;Schon et al., 2012).
Figures 7c and 7d are as Figures 7a and 7b but for the Jezero Crater simulations. The EDL conditions predicted by the model for this site are more favorable compared to those at Oxia Planum, with generally weaker winds (maximum values of~35-40 m/s between 25 and 30 km) and a smaller sensitivity to dust loading (maximum differences generally within ±15 m/s as opposed to ±30 m/s at Oxia Planum). For the density profile, similar results are obtained as at the other site, with the MCD scenario also featuring higher dust loading compared to the TES one. While at Jezero Crater the model does not predict a nighttime low-level jet, it is important to note the layer of strong winds between about 5 and 10 km in the evening hours with maximum speeds up to 30 m/s, present in both the MCD and TES simulations, that may pose a problem for the final stages of the EDL. As at Oxia Planum, the most favorable EDL time seems to be in the late morning and early afternoon. Overall, and according to the MarsWRF predictions, the EDL conditions for this mission may be less of a problem than those for ExoMars 2020.
Pla-García and Rafkin (2015) used the Mars Regional Atmospheric Modeling System (MRAMS) to predict the EDL conditions for several of the Mars 2020 landing sites including Jezero Crater during early and midafternoon. They have used a hybrid dust scenario obtained by combining the column dust opacity of the zonally averaged TES retrievals with a vertical distribution given by the Conrath-ν profile (Conrath, 1975) that is also used in MCD dust scenario. For the ranges of heights considered, MRAMS predicts maximum wind speeds of about 45 m/s, comparable to the MarsWRF forecasts (maximum wind speeds in the afternoon are~31 m/s for the MCD and~45 m/s for the TES simulations). However, the vertical profile is slightly different with the MRAMS winds generally increasing in magnitude with height up to~35 km, whereas the MarsWRF afternoon winds strengthen from the surface to about 13-15 km with smaller changes in speed from that level until the highest considered,~34 km. These discrepancies likely result from the different experimental setups, in particular with respect to the dust vertical distribution and total column opacity, and are in line with the findings of N15.

Conclusions
Together with the launch, the most critical stage of a mission to another planet is the EDL: if it is not successful, as was the case with more than two thirds of the missions to Mars, all the proposed scientific objectives will not be fulfilled and the time and money spent on the preparation of the mission will be gone to waste. In addition to electronic/software issues, schedule constraints and human errors (e.g., Bitten et al., 2006;Sauser et al., 2009), an EDL on Mars is known to be rather challenging due to the lower density, unpredictable winds, atmospheric dust loading, and rough terrain (e.g., Braun & Manning, 2007). An accurate simulation of the EDL profiles, in particular of that of the density and horizontal wind, is therefore of the upmost importance both (i) to prevent future failures, such as the recent crash-landing of the Schiaparelli EDM in October 2016, and (ii) for the design of new EDL technologies that will pave the way for the future sustained exploration of Mars. In this paper, MarsWRF is used to simulate the EDL profiles of six past (Mars Pathfinder, MER Spirit and Opportunity, Phoenix Mars Lander, MSL, and ExoMars 2016 Schiaparelli) and two planned (ExoMars 2020 and Mars 2020) missions to the planet. The model is run in a five-nest configuration, down to a spatial resolution of about 1.46 km, centered on each landing site. Two prescribed "averaged" dust scenarios are considered to investigate the sensitivity of the atmospheric profiles to the dust loading: MCD and TES. The former aims to fit most of the thermal profiles observed during the MGS mission, while the latter is obtained from the interpolation of the TES limb opacities, combined with the total column opacity from TES nadir observations, effectively averaged over L s 110°M Y24 -L s 80°MY27. These two scenarios differ not only in the total column opacities but also in the vertical distribution of dust, which generally decays faster with height in the MCD profile. In the vertical 45 levels are considered with the model top at 120 km and the uppermost level with useful data located at about 34 km AGL. This set up is justified for the purpose of this work as the most critical phase of the EDL occurs in the lowest 30 km (Way et al., 2013). The model predicted profiles are compared to the reconstructed EDL profiles bearing in mind the errors that may arise from the uncertainties in the estimation of the atmospheric wind speed, angle of attack and drag coefficient of the spacecraft.
For the six past Mars missions the model is found to have some skill in simulating the observed EDL profiles of density, pressure, and temperature (and for Schiaparelli's EDL, that of the horizontal winds as well). In line with N15, the predicted profiles are found to be very sensitive to the total and vertical distribution of dust. Out of the two dust scenarios considered, none is found to always give the best agreement with the reconstructed EDL profiles, with both likely underestimating the dust loading for heights above~15 km (except for the Schiaparelli EDL), which is reflected in the colder temperatures and higher densities compared to the reconstructed profiles. For the Spirit EDL, these discrepancies are more significant given the nearly isothermal temperature profile below 30 km as a result of a dustier environment (Withers & Smith, 2006;N15). The attribution of the MarsWRF predicted colder temperatures and higher densities than those in the reconstructed profiles above 15 km to a clearer atmosphere in the model can be confirmed by using a dustier scenario, such as one based on the Viking Landers' 1 and 2 dust opacity observations with the global dust storms removed, as shown, for example, in Montabone et al. (2006). In addition, considering more extreme scenarios, such as those representative of regional or global dust storms, would be of interest to investigate their impact on the EDL profiles. Such simulations with different dust scenarios and/or vertical grid resolutions will be left for future work. The lack of a water cycle and hence representation of the radiative effects of water ice clouds may also explain some of the discrepancies, in particular for the Pathfinder EDL. Despite not being able to simulate the finer-scale structure in the zonal and meridional wind profiles, MarsWRF successfully captures the maximum horizontal wind speed for heights between 2.5 and 9 km AMR for the Schiaparelli EDL. In particular, it shows maximum speeds of about 20 m/s in the lowest 10 km, as in the reconstructed profile. In conclusion, and despite the referred deficiencies in the MarsWRF profiles for the six past missions considered here, its predictions of the EDL profiles can be used for guidance into the EDL of future planned missions.
The model is subsequently run for the proposed landing sites of ExoMars 2020 (Oxia Planum) and Mars 2020 (Jezero Crater) missions and for the expected time of EDL (L s 19°) using the same configuration. For Oxia Planum, and for the lowest~8 km for the horizontal wind speed and~15 km for the density, the predictions obtained using the two dust scenarios are generally similar. However, above those levels they can be very different, in particular for the horizontal wind speed with maximum differences of ±30 m/s (the density ratio is within 10%). For this site the model predicts rather strong westerly winds mostly between about 15 and 35 km that, in particular in the evening and nighttime hours, can have speeds of up to 60 m/s. For the Jezero Crater, while the simulated density variations are similar, the horizontal wind speeds are generally weaker (not exceeding 40 m/s) and show a smaller sensitivity to changes in the dust vertical profile (maximum differences of ±15 m/s). However, for this site the model predicts a layer of strong westerly winds between 5 and 10 km in the evening hours with maximum speeds of around 30 m/s that may be hazardous for the EDL. For both sites, the best EDL time seems to be in the late morning and early afternoon, when the horizontal wind speed is generally the weakest with smaller shifts in wind direction.
In this paper, it is concluded that a simplified Martian mesoscale model can give a skillful representation of the observed EDL profiles for past missions and hence can be used to predict what future missions to the planet may encounter. A correct prediction of the vertical profiles of atmospheric variables such as density and winds is also important for the planned MSR missions that will involve launching a sample from Mars in order to send it back to Earth for analysis (McLennan et al., 2012;Stefaan & Beaty, 2018). This is an important step to prepare for a future human exploration of the planet. The findings of this work will be of help to those involved in the predesign phase of such missions to Mars that can be led not just by space agencies such as NASA and ESA but also by private companies or even joint private/public missions. For the design of future missions, we recommend to run a full set of simulations that thoroughly evaluate the potential dust impacts on the EDL conditions, including regional and global dust storms that will have an effect on the model predictions. This would give engineers an "envelope" of possible EDL conditions, considering both the most extreme cases that might occur, as well as some nominal and most probable cases as shown here. By doing so, the risks of each mission can be properly assessed and incorporated into the design of the mission and in the selection of the best window for entry and launch.