The Importance of Isentropic Mixing in the Formation of the Martian Polar Layered Deposits

Layers of ice and dust at the poles of Mars reflect variations in orbital parameters and atmospheric processes throughout the planet's history and may provide a key to understanding Mars's climate record. Previous research has investigated transport changes into the polar regions and found a nonlinear response to obliquity that suggests that Mars may currently be experiencing a maximum in transport across the winter poles. The thickness and composition of layers within the polar layered deposits (PLDs) are likely influenced by changes in horizontal atmospheric mixing at the poles, which is an important component of the transport of aerosols and chemical tracers. No study has yet investigated horizontal mixing alone, which may be influenced by polar vortex morphology. We show that mixing in an idealized Martian global climate model varies significantly with obliquity and dust abundance by using a diagnostic called effective diffusivity, which has been used to study the stratospheric polar vortices on Earth and to understand their role as a mixing barrier but has not been applied to Mars's polar vortices. We find that mixing in the winter southern hemisphere doubles with either an octupling of dust loading or an increase in obliquity from 10° to 50°. We find a weaker response to changing dust loading or obliquity in the northern hemisphere. We demonstrate that horizontal mixing is an important component of transport into Mars’s polar regions, may contribute to the formation of the PLDs, and presents effective diffusivity as a useful method to understand mixing in the Martian atmosphere.


Polar Layered Deposits
The polar layered deposits (PLDs) are layered ice deposits at the poles of Mars that may provide a key to understanding the Martian climate record.Understanding their formation remains an open and important question in the study of Mars's atmosphere and its evolution (Smith et al. 2020) and is currently listed as one of NASA's higher-priority goals for the scientific exploration of Mars (MEPAG 2020).The PLDs were first observed by Mariner 9 during southern early spring (Murray et al. 1972) and may represent variations in orbital parameters (Laskar et al. 2002), as well as indicate changes in transport into the polar regions throughout Mars's history.
Mars has undergone much larger changes in its orbital parameters (including obliquity and eccentricity) than Earth has, with initial predicted variations in obliquity of 15°-35° ( Ward 1973), although more recent statistical calculations suggest variations in obliquity of 5°-60°and in eccentricity of 0-0.12 over the past 10 Ma (Laskar et al. 2002(Laskar et al. , 2004)).Previous studies have linked changing orbital parameters with the layering observed in Mars's polar regions (e.g., Perron & Huybers 2009;Hvidberg et al. 2012), but attributing these changes to different processes on Mars still needs to be examined.In particular, no studies have considered the isentropic transport of dust and other atmospheric tracers in Mars's atmosphere, which is typically associated with eddy activity.This dynamical "mixing" from chaotic advection stretches and deforms material contours of tracers, along which sub-grid-scale molecular diffusion occurs.In this study, we aim to explore what role atmospheric mixing may have played throughout Mars's history, its dependence on orbital forcing and global dust abundance, and how it may have contributed to the layering we observe at the poles of Mars.
The bulk composition of the northern PLD (NPLD) is ∼95% water ice (Grima et al. 2009), but individual layers are thought to contain varying amounts of dust and CO 2 ice (see Figure 1).The southern PLD (SPLD) is also thought to be primarily composed of water ice (with 0%-10% dust; Plaut et al. 2007), although these layers sit on a partially buried ≈1 km bed of CO 2 ice (Smith et al. 2020, and references therein).There are significant complexities in understanding the composition of individual layers of the PLD, as the reflected intensity of individual layers depends on local slope, surface roughness, and layer dust and water/ice proportions, as well as illumination due to the season of observation (Laskar et al. 2002;Perron & Huybers 2009).The extent to which layering in the PLD is due to enhanced sublimation of water ice during some periods (leaving behind layers of higher dust concentration) or due to increased water ice deposition relative to dust deposition is also unclear, although recent studies hypothesize that it is likely a combination of the two effects (Courville et al. 2021).

Atmospheric Circulation
The winter atmosphere of Mars is characterized by a strong polar jet and annular potential vorticity (PV) structure, i.e., with a minimum in PV over the pole and a region of higher PV surrounding it.The annular PV structure contrasts with the morphology of Earth's stratospheric polar vortices, which are monopolar.PV is often used in the study of polar vortices and has been investigated on Earth, Mars, and other planetary bodies (e.g., Mitchell et al. 2015;Waugh et al. 2016;Sharkey et al. 2021).Mars's annular PV structure is thought to be influenced by the release of latent heat during CO 2 condensation (e.g., Seviour et al. 2017;Toigo et al. 2017;Ball et al. 2021), along with Hadley cell downwelling (Scott et al. 2020).The atmosphere of Mars experiences particularly large interannual variability, due primarily to the variability in the occurrence and size of dust storms; during particularly dusty years, the Hadley circulation may expand poleward and the PV structure may collapse to a monopole (e.g., Ball et al. 2021;Streeter et al. 2021).Previous studies have suggested that the annular polar vortex may act as a partial barrier to transport across the polar regions (Waugh et al. 2019;Toigo et al. 2020).Furthermore, modeling studies suggest that the amount of CO 2 condensing may have varied significantly with obliquity or peak dust loading, possibly with significant implications for dust and ice sedimentation over the poles (Toigo & Waugh 2022).
Previous studies have simulated the atmosphere of Mars at various obliquities and/or eccentricities (e.g., Haberle et al. 2003;Newman et al. 2005;Toigo et al. 2020).Haberle et al. (2003) showed that obliquity has the stronger influence of the two parameters on the climate system owing to its larger effect on the latitudinal variation of solar insolation.Similarly, Mischna (2003) showed that the location and stability of surface ice are most closely linked to obliquity in Mars global climate model (MGCM) simulations varying obliquity, eccentricity, and longitude of perihelion.Newman et al. (2005) showed that at low obliquity there are permanent CO 2 ice caps at both poles, lowering global mean surface pressure.At high obliquities, they found that there was a stronger and broader meridional circulation and a larger seasonal CO 2 ice cap and that this stronger circulation increased dust lifting by wind stress, providing a positive dynamical feedback.However, no studies have yet considered the impact of a finite surface dust supply in varying-obliquity simulations, which may limit the positive feedback effect as was shown to be the case in simulations of the present-day Martian atmosphere (Newman & Richardson 2015).
Several modeling studies have investigated timescales for transport in the Martian atmosphere, including Barnes et al. (1996), who looked at "ventilation timescales" and found that there was faster ventilation of the atmosphere above 1 hPa (i.e., more vigorous transport) during the dusty solstice season compared to equinox.An idealized "age-of-air" tracer initialized in an MGCM was found to be oldest in the polar lower atmosphere, and a transport barrier intrinsically linked to polar vortex morphology was identified, with the presence of smaller-scale vortices (or "patches" of PV) being linked to a decrease in "age of air" within the annulus of PV (Waugh et al. 2019).Toigo et al. (2020) used this same tracer in a parameter sweep over obliquity to identify a turning point in transport into the polar regions, with the most efficient transport occurring at approximately Mars's current obliquity across the simulations performed.
Other work has studied the impact of dust loading on the Martian atmospheric circulation.Elevated dust loading can influence features of the northern winter polar vortex by, for example, initiating a "transient polar warming" (in which polar temperatures may rise 2-10 K during a typical year or more during a global dust storm (GDS) at high altitudes over the pole and by strengthening the polar jet and steepening the poleward tilt of the jet core at low altitudes (Guzewich et al. 2016).This warming occurs via a strengthening of the overturning circulation with higher dust loading during the dusty northern winter, which warms the polar upper atmosphere through stronger adiabatic heating from the descending branch of the Hadley circulation.Toigo et al. (2020) also suggested that larger polar warmings (in which temperatures can rise up to 70 K between 1 and 0.1 Pa) occur each simulated winter for obliquities above 35°under typical dust conditions, while they are reliant on the presence of higher dust loading at present-day obliquity.An "opposite season" dust event (with peak dust loading prescribed during the southern winter, rather than during northern winter, when it is always observed) was shown to lead to a more vigorous overturning circulation compared to standard dust loading at that time (Guzewich et al. 2016).Hemispheric differences in the response to obliquity have also been noted; Toigo et al. (2020) found that "age of air" is consistently older at the north pole in all simulations regardless of the timing of peak dust loading and proposed that dust deposition would vary on half the period of a full obliquity cycle between 5°and 55°(due to an increase in transport up to obliquity ε = 25°and a decrease above this).
In this paper we examine the influence of polar vortex morphology on transport into the Martian polar regions under different orbital forcings or global dust abundances.We vary obliquity, eccentricity, and dust scaling and use a diagnostic called "effective diffusivity" to understand horizontal mixing in the polar regions in an idealized MGCM.Effective diffusivity is a geometric description of the stirring in a flow and has been well used to study mixing in Earth's stratospheric polar vortices (e.g., Haynes & Shuckburgh 2000a, 2000b;Abalos et al. 2016), as well as across hurricane eyewalls, which have an annular vorticity structure similar to Mars's current polar vortices (Hendricks & Schubert 2009;Mitchell et al. 2021).Hendricks & Schubert (2009) find two regions of high effective diffusivity (indicating strong mixing) in simulations of unstable rings of vorticity-an inner mixing region within the annulus, and an outer mixing region outside of the annulus.The annulus of vorticity itself has low effective diffusivity, separates the inner and outer regions, and acts as a partial barrier region.Alongside this modeling of rings of high vorticity, Hendricks & Schubert (2009) also studied effective diffusivity in a storm with monotonically increasing vorticity.For these monopolar vortices, similarly to Earth's stratospheric polar vortices, Hendricks & Schubert (2009) found that the center of the storm acted as a partial barrier and air was not easily mixed.
The use of effective diffusivity allows us to isolate the changes in isentropic mixing (or mixing that occurs on isentropes, surfaces of constant potential temperature) from changes in large-scale transport via the Hadley circulation.This distinction is an important one, as isentropic mixing primarily refers to mixing due to eddies and acts to flatten tracer distributions on these surfaces, yet the diabatic mixing that occurs across isentropes tends to have the opposite effect (Andrews et al. 1987;Haynes & Shuckburgh 2000a).Furthermore, given the generally decreasing concentrations of dust with altitude, the deposition of dust via mixing has the potential to be larger than that via Hadley circulation, which may in general transport material into the polar regions at higher altitudes.As yet, no studies have used the effective diffusivity diagnostic to understand mixing in the Martian atmosphere.
Here we attempt to understand whether the annular polar vortex on Mars behaves in a similar way, acting as a mixing barrier between the polar regions and midlatitudes.We also consider how effective diffusivity (and hence mixing) changes under different orbital parameters or dust abundances, and thus how it may contribute to the layering observed in the polar regions.

Isca
We use the flexible modeling framework Isca (Vallis et al. 2018) to model the atmosphere of Mars.Specifically, the MGCM component of Isca (henceforth Isca-Mars) has been developed to simulate the circulation of the Martian atmosphere at varying levels of complexity (Thomson & Vallis 2019); further processes have since been added to more accurately represent the modern Martian atmosphere (Ball et al. 2021).The processes and data that are currently available in Isca-Mars are the SOCRATES radiation scheme (Manners et al. 2017), a prescribed dust profile, and latent heat release representation associated with CO 2 condensation, as well as topography from the Mars Orbiter Laser Altimeter (Smith et al. 1999).While a full model description is given in Ball et al. (2021), we include here a brief description of the dust profile in use.Initially, zonally averaged infrared absorption column dust optical depth (CDOD) from the Mars Climate Database (MCD) "standard dust scenario" (Montabone et al. 2015) normalized to the reference pressure of 610 Pa product is inputted and converted to a surface mass mixing ratio by the SOCRATES interface within the model.This then informs the surface mass mixing ratio parameter a 0 within the modified Conrath-ν profile as described in Montmessin et al. (2004) and Ball et al. (2021), from which we obtain a vertical and latitudinal profile.The MCD scenario is built by averaging climatologies of MY 24-31 CDOD (excluding the GDS events of MY 25 and MY 28), thus giving a dust profile with "typical" dust loading (Montabone et al. 2015;Ball et al. 2021).Figure 1 of Ball et al. (2021) shows the temporal evolution of the latitudinal and vertical profiles of dust used within Isca-Mars.The latent heat release representation within Isca-Mars is a simple temperature floor scheme, wherein model temperatures are prohibited from falling below the pressure-dependent frost point of CO 2 , and is detailed in Ball et al. (2021).The atmospheric mass lost during this phase change is not considered within the model, although it has been previously shown that the dynamical effect of the associated pressure changes has less impact on the structure of the northern Martian winter polar vortex than the latent heat release associated with the CO 2 condensation (Waugh et al. 2019).
The advantage of using Isca herein is its flexibility: although it is an idealized model and several processes are not represented (e.g., radiatively active ice clouds, an active dust scheme, or CO 2 ice microphysics), we are able to modify basic planetary parameters with ease and hence isolate some of the most important processes of interest.Despite the idealized nature of Isca-Mars, the inclusion of dust within the SOCRATES radiation scheme gives reasonable agreement with reanalyses (e.g., the OpenMARS data set described in Holmes et al. 2020) of the Martian atmosphere (see Figure 4 in Ball et al. 2021, which compares temperatures and winds of different Isca configurations with the OpenMARS reanalysis data set).

Simulations
To investigate the transport and mixing in the polar regions, we run a suite of simulations in Isca, varying orbital parameters and dust scaling.The basic simulation is the full version of Isca-Mars (Ball et al. 2021, with all features "switched on"), and from there obliquity, eccentricity, and dust scaling are varied.Key fixed and varying model parameters for each experiment are described in Table 1.Obliquity (ε) is varied between ε = 10°and 50°at intervals of 5°, representative of obliquity variations over the past 5 Ma and beyond (Laskar et al. 2004), and eccentricity (γ) takes the value of either γ = 0.093 (current) or γ = 0 (to isolate the effect of varying   Area-weighted polar average (75°-90°N/S) of tracer concentration at 300 K, t sols after initialization only obliquity and to remove any influence from the timing of perihelion in our results).
The current obliquity of Mars is ε = 25°.19,and its eccentricity is γ = 0.093, a significantly more eccentric orbit than that of Earth.Perihelion on modern Mars occurs at the areocentric solar longitude (L s ) of ∼251° (Haberle et al. 2003), not long before northern hemisphere (NH) winter solstice (L s = 270°). 3Newman et al. (2005) found that atmospheric circulation is strongest when the longitude of perihelion aligns with the NH winter solstice (due to hemispheric topographic differences), suggesting that modern Mars may be experiencing particularly strong atmospheric circulation.We find that our simulations at ε = 25°and γ = 0.093 give quantitatively similar results to those at ε = 25°.19 and γ = 0.093 (not shown), and so we take this to be our present-day reference simulation.
We perform γ = 0 simulations to isolate the effect of changing obliquity, and when γ = 0.093, perihelion occurs at its modern-day timing of L s ≈ 251°.We do not perform a parameter sweep over further values of eccentricity, as it would necessitate consideration of precession of perihelion.Furthermore, it has been shown that obliquity has the larger influence on atmospheric circulation, so our focus will be on this (Haberle et al. 2003;Mischna 2003).Finally, we also vary the dust mass mixing ratio by scaling the constant reference mass mixing ratio a 0 at pressure p 0 by λ = 1/2, 1, 2, 4, 8 (Equation (5) in Ball et al. 2021).This range of dust scalings includes a peak dust loading during the dusty season that is approximately equivalent to 2× the loading of the MY 28 GDS, which was found to be sufficient to cause a breakdown of the vortex in Isca-Mars in Ball et al. (2021).
We note that by using a present-day dust seasonality (via our choice of both MCD scenario and Conrath-ν profile) we impose a dusty season during northern winter for all simulations, with much lower dust loading over the southern hemisphere (SH) winter period.Given the uncertainty in historical dust opacity, and due to the topographic forcing that enhances circulation during northern winter, this simple approach is one plausible way to systematically vary dust loading and capture the dynamical effects of greater dust abundance.This constant scaling of dust opacity throughout the year is a different approach from that taken by previous work (see Guzewich et al. 2016, who vary the timing and magnitude of northern winter peak dust loading and maintain a small dust loading during southern winter in most simulations, except their "opposite season" simulation).Note that the magnitude of their "opposite season" southern winter dust loading is still smaller than the imposed southern winter dust loading in our λ = 4, 8 simulations.We do not explore the influence of the timing of peak dust loading here, which is beyond the scope of this study.
All simulations are performed at T42 spectral resolution, corresponding to a 64 × 128 (∼2°.8 × 2°.8) spatial grid, with 25 unevenly spaced vertical sigma levels reaching approximately 0.05 Pa (∼90 km), with enhanced resolution near the surface and the model top.

Effective Diffusivity
Alongside analysis of model output fields such as zonal winds and temperatures, we calculate Ertel's PV to understand the morphology of the polar vortices. 4Despite PV being considered a tracer over short timescales on Earth (e.g., Waugh 2023), one cannot say the same in the atmosphere of Mars owing to the destruction of PV over the poles by the release of latent heat by CO 2 condensation.This process warms the polar air and reduces vertical potential temperature gradients, thereby reducing polar PV: PV is no longer conserved owing to this diabatic heating, and so it should not be considered as a tracer over the winter poles of Mars.To avoid this complication, we initialize a passive massless tracer c that is advected by the flow.The tracer is initialized with a sinusoidal profile and is allowed to freely evolve for 30 sols before being reinitialized with the same profile.Initial tracer distribution is given by Using the passive tracer, we turn to a quantity called effective diffusivity (κ eff ) to investigate isentropic mixing.This is a geometric method of measuring the mixing occurring in a flow that looks at the lengths of contours of a conserved tracer.
Following Nakamura (1996Nakamura ( , 2008)), for a monotonic tracer field c * , one can transform the tracer field into a contour-based coordinate.Suppose that we have a contour of constant tracer c * = c.The area A(c, t) occupied by the tracer c * c is such that and it evolves according to where κ is the molecular diffusion coefficient.Given the monotonicity of the tracer, it is possible to invert the function A (c, t) to c(A, t).With some intermediate steps (following Nakamura 2008), we arrive at a relationship for the evolution of the contour c, Here the subscripted angular bracket • c á ñ is the average of a field variable (denoted here by •) on the tracer contour c * = c, such  4, On a sphere, one can then consider the area A as a function of equivalent latitude f e , which is defined by is constant on the contour, note that L e reduces to the actual perimeter length of the contour.Effective diffusivity, * eff k , is defined by where L min the minimum possible length of a contour enclosing A, as given in Nakamura (1996Nakamura ( , 2008) ) and Haynes & Shuckburgh (2000a).On a sphere of radius r, L min = r 2 cos e p f.The key idea of this calculation is to transform the advection-diffusion equation for the evolution of a tracer into a diffusion-only equation.
Mixing is due to the combined effect of differential advection and turbulent diffusion-advection stretches and deforms material lines while diffusion accomplishes true irreversible mixing.Effective diffusivity combines irreversible diffusion effects with reversible advection effects, via the equivalent length (which quantifies both the deformation of material lines by advection and the interface over which diffusion can act; Hendricks & Schubert 2009).Equivalent length absorbs the effects of advection since the tracer contour is used as the coordinate (the RHS of Equation (4) demonstrates that c only evolves with contour length via L e ).If tracer contours are well stirred (i.e., contours are well stretched and deformed by advection), one sees  L L e min , and hence effective diffusivity will be large (Qian et al. 2019), as there is more interface for diffusion to act over.Within this paper, we henceforth consider the normalized effective diffusivity, * ( ) , with the logarithm taken to more clearly illustrate large variations in effective diffusivity.
An example tracer contour is given in Figure 2 to illustrate the geometry of a flow with high effective diffusivity (left; noting the large contour length over which diffusion may occur) and low effective diffusivity (right).In all figures in which κ eff is shown, we ensure that the tracer has been allowed to evolve for at least 20 sols after initialization to minimize any influence of initial tracer distribution.
Alongside effective diffusivity, the following metrics are used to diagnose the circulation: the strength of the polar jet at 50 Pa (u max ) and its latitude ( u max f ), the strength of the overturning circulation ψ ( max y ) and the latitude of its descending branch (f HC ), and the strength of the polar vortex at 300 K (PV max ) and the latitude of the maximum in PV (f PV ; this determines the width of the annulus).We show the area-weighted latitudinal average of κ eff poleward of 40°N/S ( eff 300 k ) to illustrate changes in the strength of mixing in the midlatitudes and polar regions.This more equatorward latitude limit ensures that we capture the existence of the mixing barrier (should one exist), which may extend significantly equatorward.We also briefly discuss overall transport across the polar regions by considering how much tracer is transported into/ away from the poles in a 30-sol time period.This is simply a ratio of tracer concentration poleward of 75°N/S after 30 sols of evolution (c pole 30 ) compared to tracer concentration at initialization (c pole 0 ).All metrics and diagnostics used in the analysis of our simulations are included in Table 2.A schematic of the overturning circulation shows the latitudes and altitudes of f HC , u max f , and f PV (Figure 3).

Circulation Response to Systematically Varying Parameters
The following section illustrates the circulation response to systematically varying dust, obliquity, and eccentricity.Further illustrative figures that help to visualize the atmospheric circulation across the range of parameters are also provided in Appendices A and B of this article.
The latitude of the descending branch of the Hadley cell (f HC ) has a clear dependence on obliquity, provided that obliquity is greater than ∼30°.Above this, in both hemispheres and regardless of eccentricity value, f HC sees a shift toward the pole over the winter solstice, excursing further at higher obliquity during the peak dust season (Figures 4(a) and (c)).The relationship between f HC and dust scaling is somewhat weaker, although increasing dust scaling also acts to push f HC further poleward, particularly in the winter NH (Figure 4
The strength of the overturning circulation clearly increases with obliquity and also with dust scaling, in either hemisphere and for either value of eccentricity (Figure 5).This is due to the increased aerosol heating with higher dust loading, which has been shown to strengthen the Hadley circulation (e.g., Guzewich et al. 2016).The strength and width of the overturning circulation are predicted to increase with obliquity in axisymmetric theory and simulations (Lindzen & Hou 1988;Guendelman & Kaspi 2019).
The jet on the 50 Pa surface also moves poleward with increasing obliquity (Figures 6(a)-(d)), and this relationship is clear from the lowest obliquity (ε = 10°) to the highest (ε = 50°).There is also a shift poleward with increasing dust scaling in the SH, but there is no clear response to dust scaling in northern winter, although the jet is pushed somewhat equatorward when λ = 8, reminiscent of the circulation response to a GDS (e.g., Guzewich et al. 2016;Ball et al. 2021).
The strength of the 50 Pa jet increases with increasing dust loading in both hemispheres, although it is a proportionally much greater increase in the SH (Figure 7).In NH winter, the response of the jet also saturates between λ = 4 and λ = 8, where the jet response is somewhat weaker.The dependence on obliquity is more complex.In northern winter, there appears to be little obliquity dependence up to obliquities of ∼35°.For obliquities of 40°, there is a significant drop in the strength of the jet across most of the winter period, although this appears to recover by around L s = 320°in γ = 0.093 simulations.This drop is stronger (i.e., the jet weakens more) at higher obliquities.In the SH, there appears to be a maximum in the strength of the jet at present-day obliquity, with decreasing jet strength as ε either increases or decreases.With the combined effect of high obliquity and eccentricity (Figure 7(b)), there is a drop in jet strength over winter solstice (for ε = 50°, the jet can weaken by ∼50%)-this effect appears to be mitigated in the zero-eccentricity case, wherein the weakening is much lower proportionally.
Coincident with this solstitial weakening of the jet at high obliquity in the SH is a destruction of the annulus over the winter solstice (Figures 8(b) and (h)).We see an overall increase in annularity with increasing obliquity in the SH, aside from this solstitial monopole.In the winter NH, the polar vortex is less annular in our model (Figures 8(g) and (h)), reaching a maximum in annularity at Mars's current obliquity of ε = 25°.In both hemispheres, there is a turning point around ε = 25-30 of maximum annularity.For the highest and lowest obliquities (ε = 10°and ε = 50°), the vortex is monopolar for the majority of winter.The dependence on dust scaling is similar in both hemispheres, with a clear reduction in annularity with increased dust scaling (Figure 8).
The strength of the vortex (PV max ) depends more strongly on dust scaling in the winter SH than in the winter NH and has the opposite response.In the SH, as dust scaling increases so does the maximum PV, although this response saturates around λ = 8 (Figure 9(f)).In NH winter, maximum PV appears to be independent of dust loading for λ = 1/2-2, and the larger scalings (λ = 4, 8) weaken the vortex, particularly in early to midwinter.This is consistent with our understanding of the vortex response to GDSs, although the cause of the hemispheric difference in response is not well explored.In the obliquity simulations, there is overall a strengthening of the vortex with increasing obliquity, in both hemispheres and for both values of eccentricity simulated.When γ = 0.093, there is a destruction of PV for the highest obliquities over the winter solstice (ε = 45, 50°).When γ = 0, this effect is reduced and the weakening is not as significant, meaning that the highest obliquities still have overall the highest PV throughout the winter.

Mixing and Transport Response to Systematically Varying Parameters
In this section, we consider how mixing responds to varying planetary parameters and the associated circulation responses of Section 3.1.We first consider overall transport across the polar regions throughout the winter period.Figure 10 shows the relative change in passive tracer concentration after 30 sols  ), averaged over the polar region (75°-90°N/S).Results are shown for different initialization periods covering each hemisphere's winter and illustrate an example of how dust could be transported into the polar region as a function of different parameters.While our tracer is passive and does not actively influence the circulation as does dust, one may expect that more dust will be transported into/away from the poles as the tracer is.
In the SH, transport is strongly dependent on both obliquity and dust scaling, with ∼5% lost at ε = 10°compared to ∼50% at ε = 50°at present-day eccentricity (Figure 10(b)).An increase of dust scaling by 16 times leads to tracer transport increasing from 10% to 35% (Figure 10(d)).There is a similar dependence on obliquity in the winter NH (up to two times more transport into the northern winter polar region at ε = 50°t han at ε = 10°), although there is additionally a local minimum in transport at the winter solstice when ε = 35°(Figure 10(a)).In contrast, there is no clear relationship with dust scaling in the NH, with 25%-45% of the tracer being transported away from the polar region in all simulations.This may be linked to the winter jet latitude at 50 Pa, which is less dependent on dust scaling in NH winter in our simulations than in SH winter (Figures 6(e) and (f)).The northern winter jet also undergoes a proportionally smaller strengthening (weakening) when dust scaling increases (decreases) than the southern one, which may be limiting the transport response.Indeed, this limited response may be linked to the generally weaker circulation response to increasing dust scaling in NH winter compared to SH winter.
The higher dust loading imposed by the timing of perihelion in NH winter may be saturating the circulation, which does not undergo large changes until complete vortex breakdown is reached, while the SH winter circulation is more sensitive, as initial loading is lower.
When considering the differences between the current-and zero-eccentricity cases, we note that more tracer is transported out of the northern winter polar region in current-eccentricity simulations compared to zero-eccentricity simulations, whereas the reverse is true in the SH.This result may be due to the timing of perihelion in current-eccentricity simulations, which leads to enhanced circulation during NH winter compared to simulations in which γ = 0.In the southern winter, however, we expect that circulation will be enhanced in simulations with γ = 0 compared to γ = 0.093, as southern winter occurs during aphelion (so that when γ = 0, insolation will be greater over southern winter than when γ = 0.093).Finally, we also note that in both hemispheres transport is greatest over winter solstice and tends to be slightly weaker over fall and later winter.
To understand where these concentration changes come from, we now consider solely isentropic mixing through the calculation of effective diffusivity, κ eff .Figure 11 shows a strong mixing barrier in the northern winter polar regions at lower obliquities where κ eff is small, with much higher κ eff in midlatitudes.This general structure is similar to that found for Earth's stratospheric polar vortices (e.g., Haynes & Shuckburgh 2000a), where steep PV gradients at the edge of the vortices act as barriers to mixing.By the highest obliquity (ε = 50°, Figures 11(i) and (j)), this barrier has all but disappeared, and there is mixing throughout the winter NH.There is clear mixing at the pole at low obliquities, where air is confined by the barrier.Figure 12 shows a similar picture, with a stronger and wider mixing barrier at lower obliquity than at higher obliquity.In the SH, mixing in the midlatitudes (equatorward of the mixing barrier) tends to increase with increasing obliquity, although we do not see this response in the winter NH.Eccentricity has a much smaller influence on the mixing barrier.In the SH, we note that the mixing barrier is generally narrower and the vortex less annular (as shown by the purple line approaching 90°S) in simulations with γ = 0 than their equivalent with γ = 0.093.This appears to be the opposite response to that of the NH winter polar vortex, where γ = 0 would suggest a slightly wider annulus and stronger barrier.
Similarly, we further see a transport barrier for different values of λ, the scaling parameter for dust loading.Mixing within the polar region is strongest when dust scaling is lowest in either hemisphere, likely due to the relatively wide annulus in these simulations.As dust scaling is increased in the winter NH (Figure 13, right panel), there is no clear response in either the width of the mixing barrier or the annulus.The only clear response of the annular vortex is at a dust scaling of λ = 4, wherein the vortex becomes monopolar at all altitudes, or at λ = 8, where the annulus is confined to lower altitudes than λ = 1/2-2.This is in contrast to the southern polar vortex, where the annulus tends to shrink with increasing dust scaling (Figure 13, left panel) but remains annular for all dust scaling scenarios.We do, however, see an increase in midlatitude mixing and that the mixing barrier contracts toward the pole in the SH with increased dust scaling.
Having considered a vertical cross section of the atmospheric circulation, we now show the evolution of effective diffusivity on the 300 K surface throughout the seasonal cycle in Figure 14.Effective diffusivity and hence mixing are largest in the winter low to midlatitudes for all simulations, particularly in the NH.We note that at higher obliquity (ε = 40°-50°) there is a significant increasing in polar mixing in the SH over the winter solstice (Figures 14(d), (e), (i), and  (j)).This is accompanied by a destruction of the annulus and hence a monopolar vortex at this time, as well as a more poleward jet.
Figure 15 shows the seasonal evolution of effective diffusivity as in Figure 14 but for the dust scaling experiments.We see that there remains an effective mixing barrier in the SH throughout winter at all dust loadings simulated, as well as an annular vortex and mixing within the annulus.The polar vortex in both hemispheres becomes less annular with increased dust loading throughout the winter period, although the southern polar vortex remains persistently annular at all dust scalings.We note that mixing in the SH winter becomes stronger at higher dust loadings (until it is approximately equivalent to that in NH winter at λ = 4 and λ = 8; Figures 15(d) and (e)).Despite there being an approximately equivalent increase in mixing at the equatorward flank of the SH winter transport barrier, we do not see the same solstitial destruction of the annulus and mixing barrier as we see in the high-obliquity simulations.
Finally, we also show changes in the large-scale atmospheric circulation with our changing parameter values.We show the winter mean (60 sols averaged about the winter solstice) values of f PV , f HC , u max f , and eff 300 k for all simulations, as well as their rates of change with obliquity and dust scaling (Figure 16).We see that f HC increases with both obliquity and dust scaling in both hemispheres, although the response tends to be slightly stronger in the SH.There is an additional eccentricity signal at high obliquity, with current-eccentricity simulations having a slightly more poleward descending branch of the HC than zero-eccentricity simulations in the SH.u max f lies approximately 10°poleward of f HC in almost all simulations, similarly increasing with both obliquity and dust scaling in both hemispheres, although with very little eccentricity signal in either hemisphere.In the λ = 8 simulation, u max f sits equatorward of f HC , consistent with the equatorward shift of the jet during GDS periods.However, the response of the annular vortex is somewhat less clear.There is a small increase in f PV in the winter NH with .The percentage of tracer remaining above 75°N/S at 30 sols after initialization.The concentration, given as a percentage, is given to be c c pole pole 30 0 in the NH, where c pole t is an area-weighted polar average of tracer concentration at 300 K t sols after initialization.Corresponding L s values for initialization (L s 0 ) and 30 sols after initialization (L s 30 ) are given in the legend.As the tracer is initialized with a sinusoidal profile (Equation ( 1)), typically it is transported into the southern polar region; hence, we take the inverse c c pole pole 0 3 0 in the SH, for comparability.Lower percentages indicate greater transport, as more tracer has been transported into/away from the polar region over 30 sols.Solid lines show γ = 0.093, and dashed lines show γ = 0. dust scaling, as well as a larger increase in the SH, consistent with the stronger response in f HC and u max f in the SH.The response to changing obliquity is less clear-in NH winter, there appears to be a turning point at around ε = 30°-35°, and similarly in the SH.In the SH, this is due to the solstitial monopolar PV structure that occurs at higher obliquities (Figure 14-while the vortex is annular, the annulus tends to be wider in the high-obliquity simulations; however, this solstitial destruction reduces the mean annularity over the winter period).Overall, we find significant differences in the width of the annulus (and its response to dust scaling) in the winter SH compared to its northern counterpart.Such differences are also identified in studies of the polar vortices in reanalyses and other free-running MGCMs (e.g., Mitchell et al. 2015;Waugh et al. 2016;Streeter et al. 2021), indicating that vortex morphology bears further investigation.We now consider the mean value of κ eff poleward of 40°N /S on the 300 K level ( eff 300 k ) in all simulations (Figure 16, purple lines).These values show that overall mixing in the SH is heavily dependent on dust scaling and obliquity.Indeed, with an increase in dust scaling from λ = 1/2 to λ = 8, eff 300 k is almost tripled (moreover, this is an increase in the effective diffusivity on a log scale, indicating that the proportional increase is actually larger at a 7.4-fold increase).With an increase in obliquity from ε = 10°to ε = 50°, there is an increase of an almost equivalent magnitude.However, in the winter NH, the response to changing either parameter is less clear.Increasing dust loading weakly increases overall mixing (Figure 16(e)), which remains large for all simulations.Increasing obliquity increases mixing for lower values of obliquity (ε = 10°-25°) but has less impact for values higher than this.We note that mixing is generally lower in the winter SH than in the winter NH, although at high obliquity or dust loading the greater sensitivity in the SH means that they have reached an approximately equal strength of mixing.
These changes may be summarized by considering the rate of change of each quantity with either obliquity or dust loading (Figures 16(a) and (b)).We see clearly that there is a poleward trend in both f HC and u max f with increasing obliquity and dust scale (except in the NH dust simulations) and that the rate of change with dust loading is much greater in the SH (Figure 16(b), dark blue).We also see an increase in mixing poleward of 40°N/S in both experiments, the response being stronger in the SH, although this hemispheric difference is more significant in the dust-scale experiments.However, we note that while f PV tends to move poleward in the dust-scale simulations, its rate of change with obliquity is not significant.Finally, Figure 17 illustrates how mixing may have varied over the past 20 Myr according to obliquity variations alone.In the NH, our simulations predict that mixing has been reasonably consistent across this time period, with lower mixing during the past 5 Myr, when obliquity has been at its lowest.However, in the SH, our simulations suggest that isentropic mixing has varied significantly with obliquity, and in particular with a strong decrease again during the past 5 Myr.Note that between 5 and 20 Ma there is a cutoff in mixing for γ = 0.093 simulations.This is due to the small nonlinearity seen in mixing between ε = 25°and 30° (Figure 16(d)).Figure 17(c) suggests that the evolution of obliquity may have a more important influence in the SH than in the NH.This figure should be treated as a demonstration that mixing may have experienced oscillatory behavior driven by Mars's orbital parameters.It does not show mixing evolution with dust loading or with eccentricity (as these are harder to constrain based on our simulations and knowledge of historic dust loading) and so should be treated with caution, but it illustrates that mixing could contribute to the layering observed in the polar regions, particularly the south.

Discussion
Our simulations have shown that mixing into the polar regions of Mars has a strong dependence on obliquity and average dust loading, along with a much weaker dependence on eccentricity.This stronger response to obliquity compared to eccentricity echoes the results of previous studies using an MGCM to investigate orbital forcing of the water and CO 2 cycles, as well as global atmospheric circulation patterns (e.g., Haberle et al. 2003;Mischna 2003).The relationships we have found herein differ from those in previous transport studies (e.g., Toigo et al. 2020), likely due to model differences and the processes represented.Despite these differences, this work uses a novel technique in the study of the Martian atmosphere and successfully isolates changes in isentropic mixing from changes in transport, which have not been identified by previous studies.
In their work, Toigo et al. (2020) found a maximum in atmospheric transport into the polar regions at around ε = 25°, whereas in this paper we find that transport increases monotonically with obliquity (Figure 10(a)).This difference may be linked to our representation of either polar tracer concentration as a proxy for overall transport.Model differences may also contribute here: for example, the vertical resolution of our simulations is lower than in Toigo et al. (2020), and Isca-Mars represents latent heat release via CO 2 condensation in a simpler scheme that does not update atmospheric pressure according to the mass lost.When we  Laskar et al. (2004; data obtained at http://vo.imcce.fr/insola/earth/online/mars/mars.html), and mixing changes with obliquity in (b) the NH and (c) the SH.We have chosen to include only variations based on obliquity owing to the complication of the timing of perihelion with varying eccentricity also, although we show the time series for both γ = 0.093 (solid blue) and γ = 0.000 (dashed red).
separate isentropic mixing from transport, however, we identify a weak maximum in winter-time mixing in the NH over ε = 25°-35° (Figure 16(c)).Despite these differences in the precise dependence, both studies suggest that mixing and transport changes must be considered when discussing the formation of the PLD.We see that, over the obliquity values Mars has experienced in the past 10 Ma, transport into the polar regions may almost double in the NH or increase almost 10fold in the SH based on changes in obliquity alone.Similarly, an octupling of atmospheric dust loading from a current climatological loading may induce approximately a tripling of transport into the southern polar region, as well as an overall increase in mixing poleward of 40°S by approximately twofold.
It is thought that the layers within the PLD are formed of differing concentrations of dust and water ice.We have shown herein that there is a stronger mixing and transport dependence on dust loading in the southern winter polar region than in the northern winter polar region and that mixing into the southern winter polar region is very high when dust scaling is also high.Given the positive feedback effects noted in Newman et al. (2005), whereby dust lifting increases with obliquity, one may expect even greater mixing increases with obliquity when an active dust scheme is used, and therefore our results may be considered to be conservative.
We note that mixing dependence on eccentricity is linked to the morphology of the polar vortex.In all simulations, if there is a change in eccentricity that is concurrent with a small shift in the annularity of the polar vortex, there is an equivalent change in the strength of the barrier to mixing.For example, in the SH winter, changing eccentricity from its current value to zero reduces the annularity of the vortex, and we see a concurrent increase in mixing into the polar region (e.g., Figure 16(d), solid to dashed lines).An equivalent reduction in annularity of the northern winter polar vortex similarly sees an increase in mixing, although this is a response to increasing the eccentricity rather than reducing it (e.g., Figure 16(c), dashed to solid lines).Due to the timing of perihelion just prior to northern winter solstice, both such changes are linked to an increase in insolation and a strengthening of the overturning circulation.
The mixing dependence on obliquity appears less related to the annularity of the vortex, which has a nonlinear response to increasing obliquity in both hemispheres, while mixing does not (Figures 16(c) and (d)).In both hemispheres there is an overall increase in transport into the polar regions (Figures 10(a) and (b)) with increasing obliquity.The mixing barrier shrinks and shifts poleward (e.g., Figures 11 and 12), meaning that there are simulations with an almost monopolar vortex where there may exist a strong mixing barrier or no mixing barrier (Figures 11(a) and (i)).
We have additionally identified a broad range of obliquities (ε = 25°-35°) that act as a turning point for some of the metrics studied.For obliquities greater than this, we see poleward excursions of the edge of the Hadley cell and the jet (Figures 4  and 6), as well as significant weakening of the jet around winter solstice (Figure 7).We also see that the vortex is most annular around ε = 25°-35°, and above this there is a pronounced solstitial monopole, particularly in the SH.There is also a concurrent reduction in the strength of the vortex for these high obliquities (Figure 9).Despite the existence of this turning point, some quantities respond linearly to increasing obliquity, including the latitude of the jet at 50 Pa and the descending branch of the Hadley circulation.
In contrast to the work of Guzewich et al. (2016), we have seen that the southern winter circulation is more sensitive to atmospheric dust loading in our model.Guzewich et al. (2016) showed that temperatures in the southern winter polar vortex were relatively invariant to a change in dust loading (from a dust-free scenario to a scenario with several times the climatological average).Here we have seen a stronger response to dust scaling in the winter SH than the winter NH in the following quantities: annulus width, latitude of the edge of the Hadley cell, latitude of the jet at 50 Pa, and, most significantly, the mean effective diffusivity on the 300 K surface (Figure 16).We also see substantial changes in southern polar temperatures with dust loading (Figure C1).However, these differences could be due to the magnitude of the SH dust loading explored in our simulations, which, when scaled by λ = 4, 8, is substantially larger than that of the "opposite season" loading used in Guzewich et al. (2016).
There is no guarantee that dust optical properties have remained constant throughout Mars's history (Clow 1987).Studies have shown that altering dust optical properties in an MGCM can lead to almost equivalent circulation changes to a doubling of dust loading, although this is ascribed to an equivalent change in aerosol heating rates (Guzewich et al. 2016).Within this study, we have scaled dust loading throughout the Martian year, taking a climatological scenario of the current dust evolution and multiplying it by a constant throughout the year, to mimic a change in global dust abundance.Other studies have investigated the influence of dust loading by changing the scaling during the peak season, or by changing the timing of the peak season (e.g., Toigo & Waugh 2022).We have found differing transport and mixing responses to dust loading than in this previous study, and the prescribed modern scenario of dust loading may not be representative of historical dust distributions.Indeed, Newman et al. (2005) showed that obliquity changes lead to changes in dust lifting, which then influence the atmospheric circulation.A similar study with an active dust lifting scheme would be beneficial to understand how dust loading changes with obliquity would influence mixing.Inclusion of a dust lifting and settling scheme would also further our understanding of how polar temperature changes are connected to dust settling rates, which may influence the results.For example, the greater polar insolation at high obliquity elevates polar temperatures (Figures C1 and C2), which may slow the dust settling rates at these times.Based on our modeling, we may expect this effect to be greater in the winter NH than in the winter SH (Figures C1(a) and (b)), although an active dust tracer would be useful in quantifying any changes in settling rates.
Comparing the dependence of mixing on either dust scaling or obliquity, as well as the hemispheric differences noted therein, we note that for all obliquities there is more mixing in the NH winter than in the SH winter (Figure 14).However, for the dust scaling simulations, mixing is stronger in the NH but remains relatively invariant to dust loading, whereas in the SH mixing increases with dust loading until it is approximately equivalent to that in the NH when λ = 4 (Figures 15 and 16).
Our results show that effective diffusivity is a useful diagnostic of mixing in the Martian atmosphere, noting that similarities with structures on Earth provide a useful example for comparison.While there are limitations to our study, particularly in the use of an idealized model with prescribed dust, we believe that this methodology could be a valuable tool to answer questions concerning the formation of the PLD, just as it has been in understanding the chemistry of Earth's polar stratosphere.We intend this study to be seen primarily as a proof of concept for the use of effective diffusivity as a mixing diagnostic for the Martian atmosphere.While our results may not be entirely conclusive with regard to the contribution of horizontal mixing to the formation of the PLD owing to modeling choices made herein, our findings provide a foundation for more complex GCMs to further explore these results.

Summary
We have investigated the role that isentropic mixing may play in the formation of the PLDs using an idealized MGCM.To do this, we have applied a novel technique for the Martian atmosphere, namely the calculation of effective diffusivity, a Lagrangian measure of isentropic mixing.We find significant changes and hemispheric differences with obliquity, average dust loading, and, to a lesser extent, eccentricity.Our key findings are summarized below: 1. We find generally high effective diffusivity in the midlatitudes, with a partial mixing barrier approximately coincident with the equatorward edge of the polar vortex.This structure is similar to that found by previous studies of Earth's stratospheric polar vortices.2. We see that midlatitude mixing in the winter southern hemisphere is more sensitive to obliquity or dust loading changes than that in the winter northern hemisphere.3. Despite overall mixing in the southern hemisphere being more sensitive to changes in dust loading, we see that there is a persistent mixing barrier across all simulations performed, whereas this is not the case in the winter northern hemisphere.4. The mixing barrier in the southern hemisphere is colocated with a local PV maximum away from the pole (suggesting an annular vortex), both vanishing only at high obliquities over the winter solstice, when polar insolation is greatest.5.In both hemispheres there is an increase in mixing and transport between the polar atmosphere and midlatitudes as obliquity increases, or as dust increases (only in the southern hemisphere).
Our results demonstrate the importance of horizontal mixing in the Martian polar atmosphere and link polar vortex morphology to the efficiency of the transport and mixing barrier.Changes in mixing of a passive tracer with obliquity and dust loading may have important implications for the transport and mixing of dust in the polar regions.In future research, considering the influence of precession of perihelion and the circulation response over a range of eccentricities will be important in understanding how these factors have contributed to the formation of the PLD, to complement the initial investigation here into the two illustrative eccentricity values we have chosen.Previous research has linked changes in dust loading to changing orbital parameters (e.g., Newman et al. 2005); further investigation into mixing changes with an active dust scheme would be illuminating.Furthermore, we have once again noted significant differences in the morphology of the northern and southern winter polar vortices, so an in-depth study into the causes of these differences would be beneficial.These hemispheric differences suggest that the processes involved in the formation of the PLD may play unequal roles in the northern and southern hemispheres.We hope that the methods presented herein and the hemispheric differences identified provide a strong case to suggest that isentropic mixing is an important process that may have contributed to the formation of the PLD.

Figure 2 .
Figure 2. Constant contours of tracer c * = c on a sphere of radius r.Left: an example of tracer geometry with large effective diffusivity.Right: an example with small effective diffusivity.In each, the red contour represents a contour of tracer c * = c enclosing the same area A. The blue contour is the shortest possible contour enclosing A, sits at latitude f e , and has length L r 2 cos e min p f = .In the left panel  L L e

Figure 3 .
Figure 3. Top: northern winter mean meridional streamfunction (shading; 10 8 kg s −1 ), zonal-mean zonal wind (thick contours; m s −1 ); bottom: zonal-mean PV (orange; MPVU) and effective diffusivity (blue) on the 300 K level.The latitudes of our metrics (f HC , umax f , f PV , eff300 k ) are indicated either by crosses or, if the calculation is based on an integral or mean, by bars across the region in question.The data used are from the present-day simulation (γ = 0.093, ε = 25°).

Figure 4 .
Figure 4. Winter anomaly of the latitude of the descending branch of the Hadley circulation (f HC ) for (a, b) γ = 0.093; (c, d) γ = 0; and (e, f) dust-scale simulations.Panels (g) and (h) show f HC for ε = 25°and λ = 1, and panels (a)-(f) show the anomaly from this value in the corresponding suite of simulations.The left column shows the winter NH, and the right column shows the winter SH.

Figure 5 .
Figure 5. Winter anomaly of the strength of the Hadley circulation ( max y ) for (a, b) γ = 0.093; (c, d) γ = 0; and (e, f) dust-scale simulations.Panels (g) and (h) show f HC for ε = 25°and λ = 1, and panels (a)-(f) show the anomaly from this value in the corresponding suite of simulations.The left column shows the winter NH, and the right column shows the winter SH.

Figure 6 .
Figure 6.Winter anomaly of the latitude of the jet on the 50 Pa surface ( umax f ) for (a, b) γ = 0.093; (c, d) γ = 0; and (e, f) dust-scale simulations.Panels (g) and (h) show f HC for ε = 25°and λ = 1, and panels (a)-(f) show the anomaly from this value in the corresponding suite of simulations.The left column shows the winter NH, and the right column shows the winter SH.

Figure 7 .
Figure 7. Winter anomaly of the strength of the jet on the 50 Pa surface (u max ) for (a, b) γ = 0.093; (c, d) γ = 0; and (e, f) dust-scale simulations.Panels (g) and (h) show f HC for ε = 25°and λ = 1, and panels (a)-(f) show the anomaly from this value in the corresponding suite of simulations.The left column shows the winter NH, and the right column shows the winter SH.

Figure 8 .
Figure 8. Winter anomaly of the latitude of the maximum PV value on the 300 K level (f PV ) for (a, b) γ = 0.093; (c, d) γ = 0; and (e, f) dust-scale simulations.Panels (g) and (h) show f HC for ε = 25°and λ = 1, and panels (a)-(f) show the anomaly from this value in the corresponding suite of simulations.The left column shows the winter NH, and the right column shows the winter SH.

Figure 9 .
Figure 9. Winter anomaly of the strength of the PV maximum (PV max ) for (a, b) γ = 0.093; (c, d) γ = 0; and (e, f) dust-scale simulations.Panels (g) and (h) show f HC for ε = 25°and λ = 1, and panels (a)-(f) show the anomaly from this value in the corresponding suite of simulations.The left column shows the winter NH, and the right column shows the winter SH.

Figure 10
Figure 10.The percentage of tracer remaining above 75°N/S at 30 sols after initialization.The concentration, given as a percentage, is given to be c c

Figure 12 .Figure 13 .
Figure 12.SH winter effective diffusivity for varying-obliquity simulations.Shading and contours are the same as in Figure11, but for the southern hemisphere (centered around L s = 90°).

Figure 14 .
Figure 14.Seasonal changes in mixing in varying-obliquity simulations.Shown is the evolution of effective diffusivity on the 300 K level (shading), zonal winds (black contours), and f PV (purple line).A 30-sol smoothing is applied to the model output for clarity.

Figure 15 .
Figure 15.Seasonal changes in mixing in dust-scale simulations.Shown is the evolution of effective diffusivity on the 300 K level (shading), zonal winds (black contours), and f PV (purple line).A 30-sol smoothing is applied to the model output for clarity.

Figure 16 .
Figure 16.The change in various diagnostics of the circulation with obliquity or dust scale.(a, b) Rates of change of the winter mean (60 sols averaged about L s = 270 (L s = 90) in the northern (southern) hemisphere winter) values of f HC , umax f , f PV , and eff300 k with obliquity (panel (a)) and dust scale (panel (b)) for both the northern (light blue) and southern (dark blue) hemispheres.Descriptions of the quantities shown are given in Table 2, and error bars show the standard deviation of these rates of change.(c-f) Winter average values of f HC (blue), umax f (green), f PV (yellow), and eff300 k (purple).Panels (a)-(b) show the gradient of the lines in panels (c)-(f).In panels (c)-(d), solid lines show current-eccentricity (γ = 0.093) simulations, and dashed lines show zero-eccentricity (γ = 0) simulations.Note that eff300 k is plotted on the right-hand axis of all panels.

Figure 17 .
Figure 17.A time series of winter-time mixing ( eff300 k ) according to historical integrations of obliquity.Winter mean is 60 sols averaged about L s = 270°(L s = 90°) in the winter northern (southern) hemisphere.(a) Obliquity fromLaskar et al. (2004; data obtained at http://vo.imcce.fr/insola/earth/online/mars/mars.html), and mixing changes with obliquity in (b) the NH and (c) the SH.We have chosen to include only variations based on obliquity owing to the complication of the timing of perihelion with varying eccentricity also, although we show the time series for both γ = 0.093 (solid blue) and γ = 0.000 (dashed red).

Figure A1 .
Figure A1.North polar stereographic map of PV (shading) and zonal winds (black contours) on the 300 K surface in current-eccentricity simulations, averaged over 30 sols either side of the winter solstice.ε = 25 shows absolute PV and zonal wind, and all other panels show the PV and wind anomaly from ε = 25 for the given value of ε.The purple contour indicates the latitude of maximum PV (f PV ).

Figure A2 .
Figure A2.North polar stereographic map of PV (shading) and zonal winds (black contours) on the 300 K surface in zero-eccentricity simulations, averaged over 30 sols either side of the winter solstice.ε = 25 shows absolute PV and zonal wind, and all other panels show the PV and wind anomaly from ε = 25 for the given value of ε.The purple contour indicates the latitude of maximum PV (f PV ).

Figure A4 .
Figure A4.South polar stereographic map of PV (shading) and zonal winds (black contours) on the 300 K surface in zero-eccentricity simulations, averaged over 30 sols either side of the winter solstice.ε = 25 shows absolute PV and zonal wind, and all other panels show the PV and wind anomaly from ε = 25 for the given value of ε.The purple contour indicates the latitude of maximum PV (f PV ).

Figure A5 .
Figure A5.North polar stereographic map of PV (shading) and zonal winds (black contours) on the 300 K surface in dust simulations, averaged over 30 sols either side of the winter solstice.λ = 1 shows absolute PV and zonal wind, and all other panels show the PV and wind anomaly from λ = 1 for the given value of λ.The purple contour indicates the latitude of maximum PV (f PV ).

Figure B1 .
Figure B1.Northern winter mean meridional streamfunction averaged over 30 sols either side of the winter solstice for γ = 0.093.ε = 25 shows absolute value of ψ, and all other panels show the anomaly from ε = 25 for the given value of ε.

Figure B2 .
Figure B2.Southern winter mean meridional streamfunction averaged over 30 sols either side of the winter solstice for γ = 0.093.ε = 25 shows absolute value of ψ, and all other panels show the anomaly from ε = 25 for the given value of ε.
Latitude of the descending branch of the Hadley circulation; latitude of the 10% crossing poleward of max y f PV (°N/S)Latitude of the PV maximum; latitude of PV max eff300 k Area-weighted latitudinal average κ eff poleward of 40°N/S at 300 K c pole t