Martian great dust storms: Interpretive axially symmetric models

ICARUS

55,

302-331 (1983)

Martian

Great Dust Storms:

Interpretive

Axially

Symmetric

Models

EDWIN K. SCHNEIDER Center for Earth and Planetary Physics, Harvard

University, Cambridge,

Massachusetts

02138

ReceivedOctober 20, 1982;revised March 31, 1983 The simplified theory of steady, nearly inviscid, thermally forced axially symmetric atmospheric motions developed by Schneider (1977) is applied to the study of the problem of the Martian great dust storms. A highly idealized calculation of the atmospheric response to heating concentrated in a small latitude band is carried out. Qualitatively different local and global response regimes are identified. As the heating is increased from zero, some critical value is reached at which the response jumps from local to global. It is suggested that this transition from local to global response may be related to the observed explosive growth of great dust storms. Results from the idealized model indicate that subtropical latitudes are favored for the initiation of a dust raising global dust storm, as the meridional scale of the response to a heat source of fixed intensity is largest for the heat source located close to the equator, but the surface stress in the zonal direction produced by the response increases as the heat source is moved towards the poles. Also, the steady axially symmetric Martian response to solar forcing is examined. Modification to the solar forced response due to an added latitudinally localized heat source is briefly discussed, and it is indicated that similar transition behavior to that obtained in the more idealized model is to be expected in this case also. Implications of the dynamical model for the dependence of the occurrence of great dust storms on orbital parameters are remarked on.

1. INTRODUCTION

One of the most spectacular planetary meteorological phenomena that has been observed is the Martian great dust storm. These storms appear to occur during every Martian year in the southern hemisphere spring/summer, when Mars is close to perihelion. An initial large but local dust storm seems to grow slowly in place for several days, and then rapid obscuration of a large proportion of the planetary surface occurs (Martin, 1974). The initial storm occurs in the southern hemisphere subtropics (Leovy, 1979). Martin’s study of groundbased observations indicates that after the initial storm is seen, other persistent local storms develop independently in approximately the same latitude band, and then the global phase of the conglomerate storm is entered. Global suspension of dust in the atmosphere produces significant changes in the zonal mean temperatures and thermal winds (Martin and Kieffer, 1979; Conrath, 1981). A pronounced upper level tempera-

ture maximum develops in northern hemisphere midlatitudes in dusty conditions, and application of the thermal wind relationship yields a northern hemisphere jet stream stronger and closer to the pole than that expected in clear conditions (Leovy and Mintz, 1969). The north polar region atmosphere is also much warmer in dusty than clear conditions. The global storm begins to decay soon after the widespread obscuration occurs (Pollack et al., 1979), with prestorm dustiness restored after about 60 days. Models of the onset of great dust storms are concerned with identification of the mechanisms for dust raising and dust spreading. The time evolution of a model great dust storm including a representation of dust raising has not to date been simulated .

An important input necessary for the quantitative modeling of great dust storm behavior is the threshold friction velocity, utc, necessary for dust raising at the Martian surface. u*, depends on atmospheric 302

0019-1035/83 $3.00 Copyright All rights

Q 1983 by Academic of reproduction

Press,

in any form

Inc. reserved.

MARTIAN

GREAT DUST STORMS MODELS

density, grain density, grain size, and grain cohesion, as well as the atmospheric molecular viscosity, and the gravitational acceleration. Due to the tenuousness of the Martian atmosphere, uec for Mars should be much larger than the threshold friction velocity for Earth, which is 15 cm set-* in dry conditions. Gierasch and Goody (1973) estimate uIc for typical Martian surface conditions (CO2 at 200°K and 6 mb) as 180 cm set-l. Other estimates for these conditions have ranged from 1 to 4 m set-*. Sagan and Bagnold (1975) point out that the large Martian topographic variations and the large diurnal temperature excursions are likely to produce significant variations in u*, with locations and time of day, through the dependence on atmospheric density. The friction velocity is related to the winds above the friction layer by a drag coefficient, Ca. Cd depends at least on the atmospheric static stability and on the surface roughness. Gierasch and Goody (1973) choose a value for Cd of 0.06, so that the surface wind speed corresponding to Use = 1.8 m set-i is 30 m set-‘. The correct Cd for Mars may differ greatly from typical Earth values, due to the high density of boulders that appear from Viking Lander imagery (Binder et al., 1977; Mutch er al., 1977) to be strewn about the surface. The large uncertainty in the values of uec and Cd appropriate to Mars makes the correctness of suggested great dust storms mechanisms difficult to evaluate. Gierasch and Goody (1973) examined the conditions required for the onset and persistence of a “dusty hurricane,” a local long-lived dust-raising storm driven by heating contrasts between a dusty core and a clear surrounding region. The model required an initial seed dust storm of a critical size, and treated the region surrounding the core as a balanced vortex with frictionally induced inflow in a surface boundary layer. The large-scale spreading of the dust injected into the atmosphere was not modeled. The vortex structure hypothesized by this model was not observed by the Viking

303

orbiters in the localized storms which appear in the initial phase of a great dust storm. Dust raising for the global storm is local in this mechanism. Leovy er al. (1973) proposed that the global dust storms occur due to a mechanism which involves a global mean meridional circulation, the thermal tides, and the seasonal mass exchange of the CO2 polar caps. Local dust storms driven by summer sublimation from the south polar cap are supposed to build up the atmospheric dustiness. Increased dustiness enhances the mean meridional circulation and the thermal tides by increasing the amount of the incoming solar radiation absorbed in the atmosphere. When a critical amount of dust is present in the atmosphere, the combined mean meridional/tidal circulation is of sufficient strength to raise dust over a wide area of the tropics, leading to a global storm. The large amount of dust raised by the global storm increases the lapse rate, slowing down the meridional circulation and reducing the surface winds below the saltation threshold, thus leading to the decay of the dust storm. The calculations of Leovy et al. (1973) concerning the role of the mean meridional circulation were diagnostic. The seasonal polar cap sublimation outflow plays an important role in this proposed scenario for great dust storms, confining global storm occurrence to the solstice seasons. This mechanism, or at least the proposed important role of the mean meridional circulation, is attractive since the mean meridional circulation including the contribution from atmosphere/polar cap mass exchange has presumably large seasonal variations, as indicated by GCM calculations (Leovy and Mintz, 1969; Pollack et al., 1981), and the occurrence of great dust storms is strongly dependent on season. Haberle et al. (1982) examined the response of a hypothetical axially symmetric Martian atmosphere to dust loading, and showed how dust injected in a specified subtropical latitude band would be spread

304

EDWIN K. SCHNEIDER

by the meridional circulation. The modeling approach followed the suggestions of Schneider (1977) concerning the elimination of the dependence of the results on the arbitrary free atmosphere viscosity. Absorption of solar radiation by dust was found to increase the strength of the meridional cirand the associated surface culation stresses, apparently up to some limiting value. Very large dust optical depths did not lead to a significant reduction in the strength of the surface stresses, a key point. This result suggests that dust raising may not depend only on the mean meridional circulation, since global dust storms begin to decay soon after the global phase is entered (Pollack et al., 1979). The result may be restated as saying that an axially symmetric global dust storm will be selfsustaining on seasonal time scales, contrary to the observations. However, the radiative parameterization chosen by Haberle et al. (1982) was noted to be incorrect for large optical depths at thermal wavelengths, so that a mean meridional circulation mechanism for dust raising and storm decay cannot be considered to have been ruled out. Also, one of the calculations to be discussed below indicates that the results may be different if the upper lid on the numerical model were to be raised. The dust-spreading calculations showed that dust/solar absorption strongly intensified the upper level meridional velocities and led to rapid dust coverage of the extrapolar regions, after the dust was raided to high levels. Other authors have considered the role of topography in the generation of great dust storms. Pollack et al. (1979) suggested that the thermally forced topographic winds play an important role in the growth and decay of global dust storms. Zurek (1976) deduced that topography enhanced the diurnal tidal winds in subtropical uplands over lowlands, suggesting that subtropical uplands should then be preferred locations for dust storm genesis (although this effect may be partially compensated by the dependence of the threshold friction velocity

on density). This article employs the theory of the steady, nearly inviscid axially symmetric atmospheric response to thermal forcing developed by Schneider (1977)) and examined by Held and Hou (1980) and Hou (1981) to study the role of the mean meridional circulation in great dust storm genesis. First, a highly idealized calculation of the response to a heat source concentrated in a small latitude band is performed. The heat source represents dust-enhanced absorption of solar radiation in a local dust storm. The interesting result is found that for a small heat source, the response is local, confined to the same side of the equator as the heat source, while for a large heat source, the response is global, with the resulting zonal winds and temperatures symmetric about the equator. A critical magnitude of the heat source, which depends on the latitude and depth of the heat source and, as well as the planetary radius, rotation rate, and gravitational constant is found. As the heating is increased slowly from zero, the response must discontinuously change from local to global when the critical magnitude of the heat source is exceeded. This discontinuous jump from local to global response resembles the behavior of observed great dust storms. In the context of this highly idealized model, local dust storms in the Martian polar regions produce only the local response (a global response is achieved only when the polar centered local storm extends from the pole to 52” latitude). A set of calculations of the response to a heat source of constant strength, but variable position is described. The maximum friction velocity for each case is found. A strong minimum in the maximum friction velocity occurs for heating at the equator. The largest surface stress produced by the response to the heating occurs polewards of 20” latitude. The larger the surface stresses, the more likely the threshold friction velocity for dust raising is to be exceeded. Thus, it is suggested that an axially symmetric dust

MARTIAN

GREAT DUST STORMS MODELS

raising and spreading dust storm would preferentially be generated from initial local dust storms neither too close to the equator nor too close to the pole. Some slightly more realistic calculations, in which the thermal forcing is due to solstice and equinox solar radiation, are also performed. The mass exchange between the polar caps and the atmosphere is not directly included, but acts as a heat source by keeping the atmospheric temperature from falling below the CO2 condensation temperature. The performance of the simplified theory is evaluated by comparison with the GCM results of Leovy and Mintz (1969) and Pollack et al. (1981), and with the two-dimensional simulations of Haberle et al. (1982). The simplified model predicts the gross seasonal variation of the zonal mean winds, temperatures, and meridional circulation well, indicating that neglect of the condensation flow is not likely to invalidate the model results. By far the largest magnitude of the surface stress is found for southern summer solstice forcing, suggesting that the threshold friction velocity is most likely to be exceeded in this season. Calculations are also made for the seasonal variation of a Martian general circulation model with a more realistic amount of ambient atmospheric dust than assumed in the existing GCM calculations. The model indicates that a zonally symmetric model with solstice forcing requires a very deep meridional circulation if the effects of internal friction are to be minimized. The combined problem of the response to large-scale solar forcing with enhanced dust/solar absorption in a small latitude band is then described. Potential types of transition from local to global responses in this model are identified. 2. THE MODEL

The object of the construction described below is to approximately determine, in as simple a manner as possible, the axially symmetric response to externally specified thermal forcing. It will turn out that this

305

problem reduces to the determination of the meridional scale of the response. Meridional scale in this context means the positions of internal boundaries between regions where a meridional circulation exists and regions of zero meridional circulation, rather than the more conventional concept of scale as an e-folding length. The derivation of the simplified model is based on the “nearly inviscid” assumption. There is little justification for internal friction in an axially symmetric atmosphere, while friction is expected to be important near a solid boundary. Frictional effects due to symmetric instability are, in principle, resolvable. Thus, frictional effects will be taken to be as small as possible in the free atmosphere. When there are no torques in the zonal direction, absolute angular momentum is conserved following a fluid ring. Assuming that the meridional circulation has a well-defined vertical scale, angular momentum is taken to be conserved along the upper streamline. Given that the surface friction is sufficiently strong such that the surface zonal winds are much less than the angular momentum conserving upper level zonal winds, the thermal wind equation determines the meridiomean nal gradient of the vertical temperature. By relating the vertical mean temperature to the diabatic heating via a Newtonian cooling law, it is shown that the axially symmetric response to thermal forcing has a definite meridional scale which depends on the structure and strength of the forcing. This approach was described by Schneider (1977) for an equatorial /3 plane. Extension to spherical geometry and evaluation of the validity of the approximations was performed by Held and Hou (1980) and Hou (1981). Held and Hou showed that the approximate solution compares well with numerical solution of the governing equations for a special case. The result that the structure of the response depends on the strength as well as the structure of the forcing differs qualitatively from the results of the calculation of

306

EDWIN K. SCHNEIDER

the axially symmetric response to thermal forcing in linearized models (e.g., Leovy, 1964; Dickinson, 1971). Consider, for example, the steady-state response of an axially symmetric model, linearized about a state of solid rotation, to an applied heat source of fixed structure localized in some latitude band, with parameterized radiative cooling linear in the temperature (Newtonian cooling or eddy conduction). Some linear parameterization such as Raleigh friction or diffusion is used to represent internal friction. Two properties of such a model are: (1) the structure of the resulting meridional circulation, zonal winds, and temperatures is independent of the magnitude of the applied heat source, while the magnitude of the response is proportional to the magnitude of the applied heat source, and (2) as the coefficient of the internal friction approaches zero, the atmospheric temperature approaches a state of radiativeconvective equilibrium, with the zonal wind in thermal wind balance with the temperature. The results of Schneider (1977), Held and Hou (1980), and Hou (1981), as well as the results to be presented here show that the above features of linearized models are qualitatively incorrect when the internal friction is sufficiently small. Structural changes in the response to heat sources of fixed structure but variable magnitude are possible only when the governing equations are nonlinear. In particular, in regions of an axially symmetric model in which the frictional effects are small, angular momentum is conserved following fluid rings. Angular momentum conservation is mathematically nonlinear and introduces the possibility of variable structure of the response. Nonlinearity in the heat equation does not appear to lead to structural changes in the response. Schneider (1977) has described why radiative equilibrium/ thermal wind balance is a discontinuous limit as internal friction approaches zero. The equations of motion for steady, hydrostatic, balanced axially symmetric motion on a rotating sphere are:

a4 T ag=q-g

(1)

-1 IF+_) a

ay

+ e’JIHO $ (@flow)

2Ryu + u’ am

(2)

ay’

(3)

y i =

1

= 0,

- -34

a

a(dF3

au

ay + wI?-11-2Ryv=&, a<

m

(4)

aT a

5y

(dT+rl To1 Qp +wac =c.

(5)

Here 5 = -Ho In (p/pa) where p is pressure and p. the mean surface pressure, and y is the sin(latitude). Ho is a scale height defined by Ho = R,ToIg, R, being a gas constant appropriate to the composition of the lower atmosphere, To a reference temperature, and g the gravitational acceleration. I = g/ C, is the adiabatic lapse rate (in height coordinates) at T = To, where C, is the specific heat of the gas at constant pressure. a is the radius of the planet and n is the planetary rotation rate. 4 = gz is the geopotential, u the zonal velocity relative to the surface, u the meridional velocity, and w = d{ldt the vertical velocity with respect to log-pressure coordinates. The frictional force in the longitudinal direction is represented by F, while the heating rate per unit mass is called Q. As well as describing the axially symmetric response to axially symmetric sources of heat and momentum, Eqs. (l)-(5) also describe the zonally averaged response to zonally asymmetric sources of heat and momentum. In the zonally averaged case the dependent variables U, u, w, T, and $ are zonal mean quantities, while F, and Q

MARTIAN

include contributions from the divergences of the eddy momentum and heat fluxes of the zonally asymmetric motions. Q also contains the zonally averaged diabatic heating, Insofar as the Newtonian cooling approximation for diabatic heating due to thermal radiation that is employed below is valid, the axially symmetric calculations of the response to specified heat sources (representing absorption of solar radiation, for example), may be thought of as giving either the axially symmetric response to an axially symmetric heat source or that part of the zonally averaged response to a zonally asymmetric heat source forced by the diabatic heating. The continuity equation, (2), implies that the meridional circulation (u, W) may be represented by a streamfunction such that there is no flow across streamlines. Employing the lower boundary condition, w = 0 at 5 = 0 (an approximation), the value of the streamfunction at 4’= 0 may be taken to be zero. Then separate meridional cells are bounded by the zero streamlines. Diabatic heating will be taken to be given by

as in Schneider (1977). T, is the temperature distribution that the atmosphere would have in the absence of zonal mean meridional and vertical motions, and 7 is a radiative time constant, which may be a function of latitude, height, and atmospheric composition. For an axially symmetric atmosphere, T, represents the state of radiative/convective equilibrium, so long as convection is independent of the meridional circulation. T, may be defined to include the effects of eddy heat flux divergences in the zonally average case. The heat equation, (5), using (6), may be manipulated to show that 44

307

GREAT DUST STORMS MODELS

= 0,

(7)

where A is any region enclosed by a streamline. It is assumed that the meridional circu-

lation has a definite vertical scale, H, such that the zero streamline occurs at 5 = H. This vertical scale is taken to be the vertical scale of the heat sources represented by T, or is given externally by other considerations (e.g., for comparison to a low vertical resolution GCM calculation H is taken as the height of the topmost level. H is taken as independent of y. A vertical averaging operator, -, is defined by

(-1 = ;

I,”( 14.

Taking A to be a region enclosed by a zero streamline, and assuming that between 5 = 0 and 5 = H the zero streamline is either vertical (demarking the boundary of separate meridional cells) or is a boundary of a region of no meridional circulation (streamfunction identically zero), (7) may be written as

where y1 and y2 are the extreme meridional coordinates of A. Cross differentiation of (1) and (3) to eliminate $ yields the thermal wind equation. Applying the vertical averaging operator to the thermal wind equation gives 2fl~(u(~,H)

-

4~8))

+ .&

(u2(y,H)

=

37

af

z -- -ay

-

u2(y,0))

-gv+.

when temperature

(10)

is continuous.

If temperature discontinuities are allowed at the boundary between the region of nonzero meridional circulation and radiative equilibrium, using the Margules (1906) relation, ZT F

aT j-To - + gH WI, - ay

wheref is the Coriolis parameter and 6~ is the zonal velocity discontinuity at the fron-

EDWIN K. SCHNEIDER

308

tal boundary, taken as the velocity in the Hadley region minus the velocity in the radiative equilibrium region. Continuous temperature will be assumed in the following for tractability although the frontal structure is expected from numerical simulation (Schneider, 1977). The main effect of friction is taken to be surface drag, which makes the zonal velocity at the ground small enough that U(Y,O) + U(Y,H).

(11)

Then from (10) and (11) the zonal wind at the top of the meridional cell is related to the vertical mean temperature by

--1 aT To a~ __-

1

1 Y(Y2 - Yo2)(2 - Y2 - Yo2) (1 _ y2)2 9 R [

where

-.gH

R= (flu)2

U(Y,H) 2RaU3-q

1

= The zonal winds at 4 = H are taken to be determined by the condition that in regions where there is a meridional circulation, 5 = H is a surface of constant angular momentum. This criterion represents the effect of vanishing friction in the free atmosphere, but nonvanishing meridional circulation. Since 5 = H is assumed to be the upper streamline of the meridional circulation, w = 0 there. Equation (4) with w = 0 and F, = 0 implies either u = 0 or 4YJJ)

=

MY2 diy--g

-

Y$)*

(13)

As u(yo,H) = 0, the angular momentum at (y. ,H) is equal to the angular momentum of the ground at that latitude. y. is taken to be the meridional coordinate of the zero streamline dividing the rising branches of adjacent meridional cells. The rising branch transmits the angular momentum of the surface to the upper branches of the meridional circulation. Inserting u(y,H) from (13) into (12) gives the meridional gradient of the vertical mean temperature in a region of nonzero meridional circulation as

(15)

If y. is known, (14) gives T(y) up to a constant of integration. To obtain closure of the problem in terms of T and T,, (9) is used. Taking (l/T) - (l/To), (9) may be written YIT -T

I

YI

2iRyU(Y,H) (1 +

(14)

Ldy=O

To

(16)

when (a) 7 = 7oe -W or (b) r = constant and aTtat = aTJag and (a/ay)(aTJag) = 0. Since geopotential must be a continuous function of pressure in the interior of the fluid, it is easily shown from (1) that T must be continuous. Then at the boundaries, y;, between regions occupied by meridional circulation and regions of no meridional circulation (where T = T,), the boundary conditions T = T, at yi

(17)

are obtained. Note that T = T, at an actual internal meridional circulation boundary is not implied by (17). A first-order discontinuity (a front) is allowed to occur at the boundaries of regions in radiative equilibrium and regions of meridional circulation, so long as the front is not vertical. In the case of the frontal interpretation, the motion along the zero streamline defining the front is discontinuous across the front. Additionally, the zonal wind at 5 = H is discontinuous at the yi and the meridional velocity is discontinuous at 5 = H in this construction. Figure 1 schematically depicts the structures of the meridional circu- :on, temperatures, and zonal winds that have been constructed above. Equations (14)-( 17) then determine the meridional scale of the steady axially sym-

MARTIAN GREAT DUST STORMS MODELS

0 -1

ys

0

yo

YN

1

Y-

FIG. 1. Schematic representation of construction described in Section 2. Cross-hatched area is region occupied by nonzero meridional circulation. As and AN are regions where t is dynamically influenced by separate meridional cells. Solid line with arrows is zero streamline. J+, and ys are meridional boundaries of region influenced by meridional circulation. Value of angular momentum at 5 = H between vN and ys is value of angular momentum at surface at yO.

metric nearly inviscid response to thermal forcing, as well as F and u(y,H). Note that the solution of (14)-(17) depends only on R and Fe,. There is no dependence of the meridional circulation scale, ?‘, or u(y,H) on the radiative time constant, 7, or the static stability, u = (aT/a[) + TT/T,,. However, r and o enter into the determination of u and W.

For simple representations of ir’,, (14)(17) may be solved analytically. Graphical solution is also a useful tool (Held and Hou, 1980). In Section 3, one of the simplest possible cases is considered, namely the response to b-function heating at some latitude. The graphical approach to describing the solution is employed. In Section 4 the uniformly dusty atmospheric response to large-scale solar forcing is discussed, using analytic and numerical solutions. The response to large-scale solar forcing with a superposed &function heat source is then considered. Implications of the solar forced results for the dependence of Martian great dust storms on orbital parameters are mentioned in Section 5. 3. RESPONSE TO A CONCENTRATED SOURCE

HEAT

Martian great dust storms are initially lo-

309

cal, with most suspended dust concentrated in local storms in a restricted latitude band (Martin, 1974), as discussed in Section 1. In terms of the radiative parameterization, (6), enhanced absorption of solar radiation due to dust has two effects relative to a clear atmosphere: the radiative equilibrium temperature , T, , is raised and the radiative time constant T, is shortened. Haberle et al. (1982) examined the axially symmetric response to a fixed dusty column of optical depth 5 in the visible that was located between 15 and 35”S, with the distribution of incoming solar radiation corresponding to southern hemisphere summer solstice. They found that dust/solar absorption in the restricted latitude band produced a significant global scale response when compared to a calculation of the response to the same incoming solar radiation, but in a dust-free atmosphere. The above-described calculations stimulated the following more comprehensive but less realistic calculations of the steady axially symmetric response to a concentrated heat source. The thermal driving for this highly idealized calculation is chosen to be given by

& _ TA To - 0 +

4WY

-

Yo),

(18)

where TA is a constant, S(y - yo) is a delta function, and q is called the strength of the heat source. T is taken to be a constant. The heat source is assumed to be more or less uniformly distributed between the ground and 5 = H. Thus, H in (15) is the depth of the heat source. The infinite elevation of T, at y0 is of course an unrealistic representation of the effects of dust/solar absorption, as is the restriction of the heating to a latitude band of zero width. However, given that dust/solar absorption raises T, by an amount AT, the width of the latitude band occupied by dust, Ay, may be related to 4 by Go

AY =

E

if r is the same constant

(19) in the clear and

310

EDWIN K. SCHNEIDER

dusty regions. The form (18) was chosen for reasons of conceptual and computational simplicity. By taking the background radiative equilibrium temperature as latitude independent, the response to the concentrated heat source is isolated. If q = 0, the atmosphere is in a state of no motion relative to the planetary surface; any deviation from the state of no motion is in response to the heat source. Modification to the methods and results of this section due to the inclusion of a more realistic latitudinally varying background T, is the subject of Section 4. The information gained by solving the highly idealized cases represented by (18) will turn out to be useful. The computational simplifications due to the choice (18) are that the boundary conditions (17) are simplified, and that the appropriate value of y. in (13) and (14) is the latitude of the applied heat source. The structure of the meridional circulation response to (18) in a statically stable atmosphere will be two thermally direct meridional cells with coincident rising branches with infinite vertical velocity at yo, as the response to the heat source must have continuous f. The latitude of the heat source, yo, is taken as 20 for purposes of discussion. The boundaries of the region occupied by the meridional circulation are called yN , where yN Z y. and ys, where ys c: yo. As the zero streamline separating the adjacent meridional cells must occur at yo, the value of the angular momentum of the surface at y. is taken to be that transmitted to the upper branches of the meridional cells, so that u(yo ,H) = 0. Substituting Fe from (18) into (16), the relation YNTTA -dy=q YS To

I

(20)

is found. Equation (14) is integrated to find

To

T -=--

1 y4 + Yo2(Yo2 - 2) R [ 2(1 - Y2)

conditions (17) %N) To

TA

To

To

(22)

Application of (22) at yN gives c

=

5 + i YN4+ Yo*(Yo* - 2) To

R

2(* - yN2)

1

(23)

and application of (22) at ys determines ys in terms of yN by YS = -YN

(24a)

or YN*+ YS*- YN2YS2 + yo2(yo2 - 2) = 0.

(24b)

Substituting (21) and (23) into (20) produces the equation relating the meridional boundaries of the response to the strength and latitude of the heat source: YN3- YS3+ YN - YS 6 2 1

YN4

+

+

Y02(Y02 -

2’)

1 - YN2 - i (1 + y0*(y0* - 2)) (In (s) - ln (Z))

- Rq = 0.

(25)

Substitution of ys from (24a) or (24b) into (25) then yields an equation for yN in terms of the parameters Rq and yo. For fixed y. the scale of the response depends only on the product Rq, essentially the vertically integrated strength of the heat source when g, R, and a are held fixed. Spurious solutions may be obtained from solving (24a,b) and (25) which would not be found from integration of (l)-(5). Two additional conditions which must be met are Y* T - TA -dy>O

I

for all y* < y. and 1+c,(21)

and C is a constant of integration that is determined by application of the boundary

TtYs) =-=-.

YS

To

T - TA ----0 To

(26)

(27)

MARTIAN GREAT DUST STORMS MODELS for some region ys < y < y. If (26) is violated, then the integral is zero for some y* < y. and the region between ys and y* is occupied by a closed meridional cell. The choice made for u(y,H) is then inappropriate in this region. If (27) is violated, then if the static stability, (T, is positive, there is ascending motion near using the approximate heat equation

311

FIG. 2. Graphical representation of equal area argument, (20), for S-function heat source at yO.See text for discussion.

(28) From continuity, the meridional motion near ys and 5 = H is then towards the heat source, and the choice of u(y,H) is again not appropriate. Conditions (27) and (28) have to be modified when the b-function response is combined with the response to large-scale variation in T, . The behavior of the solutions to the system (24a) or (24b) and (25) is difficult to deduce, as (25) is a complex algebraic/transcendental equation. In the special case y. = 0, so that only (24a) applies, and for yN < 1, a closed-form solution YN =

(%q)“’

(29)

is obtained. This is the case treated by Schneider (1977). As the strength of the heating increases in (29) the scale of the meridional response grows continuously. The speed of the jet streams (the maximum westerlies) increases and the jet streams move polewards with increasing Rq. The qualitative properties of the general solutions to (24a) or (24b) and (25) are most easily deduced from a graphical representation of the equal area argument, (20). From (21) T/To is symmetric about the equator with maxima at y = +yo and a relative minimum at y = 0. Figure 2 is a sketch of the behavior of (T - TA)lTo (the temperature curve) from (21) for some y. > 0. Changing the constant of integration in (21) corresponds to a vertical displacement of the temperature curve with respect to the y axis. Various choices of the position of the y axis relative to the temperature curve are given by the horizontal lines in Fig. 2. Fix-

ing the y axis in some position relative to the temperature curve, the zero of the temperature curve for y > y. is yN while the possible positions of ys are the zeroes of the temperature curve for y < yo. When more than one ys is allowed, the possible ys are labeled ysi with index i increasing with increasing distance from yo. The equal area condition, (20), can be used to evaluate q for a particular choice of position of the y axis. q, which must be positive for the construction to be valid, is the area enclosed by the temperature curve and the y axis between ys and yN. Area is negative in those regions enclosed by (T - Td/To < 0 and positive otherwise. Referring to Fig. 2, the behavior of the possible equilibrium solutions as the y axis is moved up and down with respect to the temperature curve will be discussed. Due to the boundary conditions, (22), the y axis must lie below the position labeled (a). Four distinct types of cases are revealed by the graphical approach. When the y axis lies below the position labeled (c) (i.e., when (F - TA)/To > 0 at y = 0) there is only one possible ys = -yN and the area enclosed by the temperature curve between ys and yN is obviously positive. For this case, the strength of the heating necessary to maintain the atmospheric circulation increases as the y axis is moved downwards from (c). These solutions imply strength of heating q > qI, where qI is the value of q at axis position (c). The zonal winds at 5 = H and the mean atmospheric temperature below 5 = H are symmetric about the equator in this case, with west-

312

EDWIN K. SCHNEIDER

erly u(y,H) for (y[ > y. and easterly u(y,H) for IyJ < yo. Temperature maxima occur at fyo, and a relative minimum of i; is found at the equator. The meridional circulation for this case consists of two cells, asymmetric about the equator, one extending from yN to y. and the other from y. to -yN. The meridional motions near [ = H are everywhere directed away from the latitude of the heat source by (27). As the strength of the heat source increases the extent of the meridional circulation increases. This type of symmetric about the equator response of T to a heat source in one hemisphere will be called global. When the y axis lies between (a) and (c) in Fig. 2, there are three equilibrium solutions possible at first sight, corresponding to the southern boundary of the meridional circulation being at ysl, ys2, or ys3, each corresponding to a different value of q. By (27), the ys2 solution is rejected, since 1; approaches TA from below as y approaches ys2. Additionally, the ys2 solution is inertially unstable, as the angular momentum at 5 = H decreases from ys2 - E to ys2 + E. The construction is self-consistent for the ys2 solution only if negative static stabilities are allowed. It turns out that the ysl solution is allowable, while the ys3 solution is sometimes allowable. The case where ys = ysl , the southern boundary of the meridional circulation being between y. and the equator, has positive area under the temperature curve by inspection, and satisfies conditions (26) and (27). The strength of the heating in this case increases as the y axis is moved, from a value q = 0 at position (a) to a value 911 = (A)ql at position (c). The response of T and u(y,H) to the heat source is local, being confined to the same hemisphere as the heat source; thus, this class of solutions will be called local. As q + 0, yN and ys approach y. . As q is increased from zero to qll , the extent of the meridional circulation increases until a maximum extent of the local response is reached, ys = 0 when q = qII. u(y,H) is easterly for y < y. and west-

erly for y > ~0. When ys = ys3, the solution is similar to the global solution described above, except that (T - TA)/To < 0 at y = 0. It is obvious that there is a position of the y axis, labeled (b) such that for the y axis between (a) and (b), q < 0, eliminating such solutions. Therefore, as the y axis moves from (b) to (c) in this case q increases from zero to qr . The constraint (26) further limits the solutions of this case, since for the axis in position (b), the area between ys3 and ysl is negative by inspection. The minimum q for allowable global solutions, qIrI, is found when the area between ys3 and ysl is zero. As ysl > 0 will be found in evaluating 4111, it is clear that qIII < qII. That is, there are certain values of q, qIII < q < qII , for which two equilibrium solutions, one local and one global are possible. Consider the response of the model at fixed R and y. as q is increased slowly from zero. Initially, the response is local, as shown above, since qIII > 0. The scale of the response increases gradually, remaining local, until the value q = qII is reached. When q = qII, ys = 0. An infinitesimal increase of q to q = qII + E must then cause a discontinuous change in the character of the response from local to global. The resulting T becomes symmetric about the equator, and the scale of the response increases dramatically. If the increase of q is interpreted as an increase of net dust/solar absorption due to area1 expansion of a local dust storm, and the scale of the response is interpreted as a measure of the dust-transporting capability of the large-scale response, then the model predicts that dust transport will remain localized to the area of the dust raising for a small storm, and that a global spreading of dust will discontinuously begin when the dusty area reaches a critical size. This critical area of the local dust storm can be estimated from (24b) with ys = 0, (25), and (19). When ys = 0, (24b) gives YN = YO m.

(30)

MARTIAN

Insertion gives

GREAT DUST STORMS MODELS

of ys = 0, into (25), using (30),

313

main local. Consequently, global dust transport due to the large-scale response to dust/solar absorption in an initially local storm is favored for the local storm in extrapolar regions. The critical local storm size does not depend explicitly on the static (31) stability or the radiative time constant; however, there is an implicit dependence as the relation for the critical strength of the through AT. heating, qc , at which the transition from loA heuristic argument based on the frontal cal to global response must occur, as a structure of the response to the heat source function of R and y. . Transition from local is useful in explaining the discontinuous beto global response will also occur for fixed q havior of the scale of the response. The reas H is increased. That is, it is the total gion in which the atmospheric temperature strength of the heat source, Hq, that is the is dynamically determined as the remote rerelevant quantity. Given qc , (19) is used to sponse to the heat source and the regions in find the width of the latitude band covered which the atmosphere is in radiative equiby the local storm at transition, and hence librium may be considered to be occupied the critical area of the local storm at transiby air masses of differing characteristics. tion is determined. Using the Martian values g = 3.72 m2 set-l, fl = 7.08 x 1O-5 These air masses are separated by fronts in set-i, a = 3.39 X lo6 m, TO = 200”K, and the upper levels, as shown in numerical assuming H = 2.5 x 104 m and AT = 50°K simulations of cases with the heat source at as in Gierasch and Goody (1973), the criti- the equator in Schneider (1977). Considering the fronts as zero-order discontinuities cal area for transition from local to global (discontinuities of temperature, but not response, A,, was determined and plotted pressure), the slope of the front is given as a function of y. in Fig. 3. Note that the approximately by the relation response to heating at the equator is always global, while there is a cutoff latitude, 52”, 62 fT6u (32) such that if the whole area from 52” to the 6y=g6T pole were dusty, the response would redue to Margules (1906), where Sz/Sy is the ratio of the vertical to the horizontal displacement of the frontal surface, and Su and ST are the changes in zonal wind and temperature taken in the same direction across the front. Since the zonal winds in the dynamically affected region are given by the angular momentum conservation argument, and the zonal winds are zero in the radiative equilibrium regions, (32) relates the slope of the front to the Coriolis parameter, I I I .l .2 .3 .4 .5 .6 7 .6 .9 1.0 f, and ST. Differences will be taken as dy7 yonamically active region minus radiative FIG. 3. The critical area of a localized dust storm for equilibrium region. Since u(y,H) is positive transition from local to global response, A,, as a funcfor jy[ > yo and negative for Iyj < yo, the tion of the sine of the latitude of the storm, yO, for sign of 6u is + at yN, and + or - at ys, Mars-like parameters. Beyond the latitude j indicated depending on whether (ysl is > or < yo. by the dashed line, A, is greater than the area between Then for the front near yN, the slope is posiY and the pole, indicating that only the local response is possible. tive or negative when ST is taken as posi-

314

EDWIN K. SCHNEIDER

tive or negative.

For either case (although case), the front near yN is Statidly stable, with warm air Overlying cold air. Similarly, a front near -yN is statically stable, (global response case), since SU is positive. When the response is local, ys is on the same side of the equator as yN and SU is negative. Equation (32) then gives a statically stable front regardless of the sign of 6T (again 6T > 0 is the relevant case). However, if ys were between the equator and -yo, 6~ would be negative, and f would be negative. Equation (32) gives a statically unstable frontal slope in this case, with cold air overlying warm air, independent of the sign of 6T. These cases are illustrated in Fig. 4. Thus, the discontinuous behavior of the response to the localized heat source can be explained as static instability of the front separating the region affected by the heat source from the unaffected region. The frontal explanation is equivalent to condition (27). 6T > 0 is the relevant

4. NUMERICAL

Numerical

APPROACH

Model

The semianalytic and graphical approach illustrated in Section 3 to finding the approximate axially symmetric nearly inviscid response to specified forcing is unwieldly for all by the simplest choices for T,. Even for the b-function forcing studied in Section 3, spurious solutions are obtained due to

-0

0

YO

FIG. 4. Schematic illustration of the frontal instabil-

ity argument given in Section 3. Sign of u(y,H) in dynamically active region is circled. Temperature denoted by W, warm, and C, cold. (a) Gl’positive. (b) 6T negative. Cold overlying warm is statistically unstable.

the integral approach, which must be argued away by consideration of the local balances in the zonal momentum and heat equations. Additionally, no information is given by the solution concerning the strength of the meridional circulation, although the meridional heat flux due to the overturning circulation is known (Held and Hou, 1980). In order to estimate the surface stresses, an important piece of information in the consideration of the generation of global dust storms, the mass flux of the meridional circulation is needed. Thus, a numerical approach which preserves the information contained in the integral approach but is capable of finding the meridional mass flux is useful. A vertically integrated baroclinic model is the simplest in the possible hierarchy of numerical models which can produce the desired flexibility. The integrated approach is developed in this section. The model is then used to examine the dust raising properties of the 6function forcing treated in Section 3, as well as the axially symmetric response of the Martian atmosphere to more realistic solar forcing. The model is developed from the approximate thermal wind equation, (12), the zonal momentum equation applied at 5 = H, the approximate vertically integrated heat equation, (28), and the mass continuity equation, (2). Assuming that the meridional velocity may be represented by

e-~“HoM({‘)d{‘)

d{.

(35)

The choice of M is arbitrary, but the numerical results of Schneider (1977), Held and Hou (1980), and Haberle ef al. (1982) indi-

MARTIAN

GREAT DUST STORMS MODELS

cate that the near-surface branch of the meridional circulation is confined to a thin boundary layer of depth 4 Ho. The vertical variation of the meridional velocity in the upper branch of the meridional circulation is taken to produce a height-independent meridional mass flux. This is an approximation of the structure found by Haberle et al. (1982). Then M = e’;~Ho (36) and _Ho*

,.,, =

eH/Ho

H

_

1 _

w

(37)

Ho

(4y,W

-1

-

2fly)

= e-H’HoFX(y,H). For axially symmetric calculations, leigh friction form is chosen for F,

(38) a Ra-

e-H’HoFx = -cu(y,H).

(39)

A nonzero damping coefficient, c, is found to be necessary for convergence of the numerical solutions with the methods employed here. F, may also include a specified forcing chosen, for example, to represent eddy momentum flux convergences. Hence, the numerical model can be used to study the zonally averaged atmospheric circulation. The vertically integrated model then consists of the system of equations (12), (28), (34), and (37)-(39). These give three equations for the dependent variables dy,H),

F’, and V,

2fiy u(y,H)

(

1 +

~Y,W

2&ZVl

- y*

= -g-g, _Ho*

H

eH/Ho

_

1 _

1I.

(12)

(V(y) vi=g

_ = v,

-1

- 24

= -cu(y,H).

_

(42)

and for finite a F/‘/ayat y = ? 1, u(y,H)

= 0 at y = + 1.

(43)

The approach of Leovy et al. (1973) is used to find the surface stress in the zonal direction. The zonal momentum equation, (4), is written in flux form and integrated vertically from 4 = 0 to { = H. Then the divergence of the net meridional flux of relative zonal momentum is balanced by the surface stress. H

A. ((1 - y*)~&‘o) ay

U-

= -u*Iu*l

(44)

gives the component of the friction velocity in the zonal direction, u,, assuming that the vertically integrated internal frictional force in the zonal direction is negligible. In order to find U, in terms of the model variables u(y,H) and V, the vertical variation of u is taken as u(uJi)

= jj

u(y,H)

(45)

and the contribution to the integral in (44) from the surface branch of the meridional circulation is neglected. Substitution of (33), (36) and (45) into (44) gives the relation for u,,

=

H 2aw

a aY

((1 - Y*)u(Y,H)VY)). (40)

(41)

The static stability, V, in (40) is latitude independent. The heat equation with a meridionally varying static stability may be derived from the flux form of the heat equation, bringing cr inside the y derivative. The appropriate boundary conditions when there is no mass exchange across 5 = 0 are

-u*Iu*l

!? d

Ho 1 a ay

MYJO

V(y) = 0 at y = ? 1

.

The contribution to m from the boundary layer return flow is neglected. Since w = 0 at 5 = H, the zonal momentum equation, (4), at 5 = H becomes V(Y) (; $

V(Y) (f 6

315

The friction velocity,

(46)

u,(y), is usually re-

EDWIN K. SCHNEIDER

316

lated to the surface geostrophic u( y ,0),by a drag coefficient Cd, u,(Y)

=

cdu(Y,o).

velocity, (47)

Haberle ef al. (1982) use Cd = 0.03 (frostfree terrain). Equation (46) will only be valid for Cd order 1, in which case the model is self-consistent. For Cd 4 1 a Selfconsistent vertically integrated model should use (10) rather than (12) and determination of the additional variable u(y,O) from the approximate form of (44), using (47). Knowledge of the depth of the surface boundary layer is needed as well as Cd for a self-consistent model when Cd 4 1. The point of view will be adopted here that while the magnitude of U, calculated from (46), Cd order 1, may be incorrect, (46) represents a first approximation to the relative variations of u,. Relative variations of U, will be interpreted as indicative of relative dust-raising capabilities. There is in reality some Iu*~[which must be exceeded for dust raising. A result from this model with larger [u,[ than another result will be interpreted as being more likely to exceed Iu*~I,thus indicating a forcing more likely to lead to dust raising. Equations (12), (40), and (41) with boundary conditions (42) and (43) are finite differenced and the two-step iterative procedure described by Schneider (1977) is used to find the solutions, given the externally specified parameters R, a, g, H, Ho, CT,T, c, and T,. The grid used has constant Ay. u(y,H) and V are defined at points which include the equator and poles, while the points at which T is calculated are offset from the velocity grid points by Ay/2. Grid spacing Ay = 0.05 is used in the results presented below. When c = 0 in (41), the existence and structure of solutions to the one layer model are given by the equal area argument of Section 2, as shown by example in Section 3. The finite difference numerical model, however, produces a noisy solution for c = t and eventually diverges as c is further decreased, where c depends on the

values of the various externally specified parameters. In the numerical calculations, a latitude-independent value of c is used, and results are shown for c = E, where 6 is determined by trial and error. The value of 6 is resolution dependent, and no argument has been developed for a priori determination of E. As the exact solutions are known for the highly inviscid nonlinear system, the results could prove useful in evaluating numerical methods. It is relevant to point out that there is a necessity for a minimum internal friction in a vertically resolved twodimensional model. Gierasch (1975) pointed out that in a steady two-dimensional calculation, there is no net latitudinally integrated vertical flux of angular momentum allowed at any level, and this condition may be used to approximately determine the amount of friction necessary for a given choice of external parameters. The same sort of approximations to the vertical structure of the zonal wind and meridional circulation as used above ((33) and (45)) may be used to estimate the meridional circulation contribution to the vertical angular momentum flux, and therefore the vertical flux of angular momentum due to small-scale processes necessary for a steady solution. This approach is successful in predicting the minimum viscosity needed in the cases treated by Schneider (1977) and Held and Hou (1980). It will be shown below that for a certain class of T,‘s and fixed R, there is also a dynamical necessity for internal friction. The calculations with the one layer model used the Martian parameter values R = 7.08 x lop5 set-i, a = 3.35 x lo6 m, g = 3.72 m* set-i, Ho = lo4 m, and To = 200°K. Response to &Function

Forcing

The response of the one-layer model is found for constant ?‘, except at one grid point. That is Tci = To + AT@(i - I)), where

(48)

MARTIAN

GREAT DUST STORMS MODELS

This is a representation of the &function forcing of Section 3 for the finite difference model. The response at fixed Z as AT is increased from zero resembles the local response at small AT and the global response at large AT, as predicted in Section 3. An example of the transition behavior of the vertically integrated model to concentrated forcing at y. = -0.375 for two different values of the Raleigh friction coefficient, c, is shown in Figs. 5 and 6. The stability coefficient is taken as (+ = 1.5 x 10e3 “K m-i and H = 11.6 km in this calculation. The condition u = 0 at y = y. is represented by replacing the zonal momentum equation by the specification resulting from angular momentum conservation, z4= 5.08 m set-l at y = -0.4 and u = -4.65 m set-i at y = -0.35. (One source of numerical instability as c + 0 is that the model loses track of the u = 0 at y = y. condition. The imposition of this condition allows convergence at a much smaller c.) The two values of c are c = 5 x lop8 set-I (Fig. 5) and c = 10m7set-i (Fig. 6). A reasonably sharp transition between the local and global behavior is found for c = 5 x lo-* set-i as AT is increased from 12 to 14°K. The transition is apparent in both the jump in the meridional velocity at the equator and the appearance of a westerly jet in the northern hemisphere. The effect of the friction is noticeable in the lack of discontinuities in the zonal winds and in the asymmetry of the zonal winds about the equator in the global solution (AT = 14°K). The jet in the opposite hemisphere from the heat source is weaker and closer to the equator than the jet in the same hemisphere as the heat source. For the larger value of the friction, c = 10m7set-I, the transition from local to global behavior is more gradual. As friction is increased, the transition from local to global behavior becomes less well defined, although the two types of limiting behavior are still found.

317

Calculations of the relative dust-raising capabilities of the large-scale circulation associated with dust/solar absorption in a localized region were performed with H = 2.5 X 104 m, 7 = 3.551 X IO5 set, u = 3 X 10e3 “K m-l, and AT = 50”K, using Fe from (48). The latitude of the heat source was varied and the maximum Iu,I for each calculation was found. The maximum Iu,[ occurred close to the latitude of the heat source except when the heat source approached the equator. The relationship found between the maximum Ju,[ and the latitude of the heat source is plotted in Fig. 7. The interesting feature of these results is the strong decrease of the maximum Iu,I as the heat source moves towards the equator. Dust initially present in an extratropical latitude band produces a large-scale response with significantly larger friction velocity than dust in a latitude band close to the equator. If the concentrated T, elevation is interpreted as due to dust/solar absorption in an initially local dust storm, the results of Section 3 (Fig. 3) and the vertically integrated model (Fig. 7) suggest that the large-scale response to this heat source favors growth to a global dust storm by dust transport and dust raising when: (1) The initial dust storm is not too close to the pole. The scale of the response favors global dust transport when the dust is initially located in the tropics or subtropics. (2) The initial dust storm is not too close to the equator. The friction velocity response increases strongly as the initial storm is moved away from the equator. An initial storm of fixed dust/solar absorption would induce friction velocities closer to the critical friction velocity needed for dust raising as the storm is moved away from the equator. More Realistic Thermal Forcing

While the analysis of the symmetric response to a pure a-function forcing described in the previous sections suggests a mechanism that divides local from global storms and gives a sharp transition from lo-

EDWIN K. SCHNEIDER

318

-‘4,,,,,,t, -.7 -6

-5

-.4 -.3 -.2 -.I

0

.l

2

3

.4

.6

.5

7

Y-

b 1.0

.5 z % \ EO >

14

zz&1

I-

-.5

I

I

1

-.6 -.5 -.4

I

I

I

1

-.3 -.2 -.l

0

.1

I

1

1

I

.2

.3

.4

5

Y-

FIG. 5. Transitional behavior of the vertically integrated numerical model to elevated T, at y0 = -0.375 for c = 5 x lo-* set-‘. (a) u(y,H). (b) VCj)e “‘“0. Curves labeled with AT, size of Fc elevation, in “K.

cal to global behavior, the results cannot be directly extrapolated to the situation that is of real interest. Although a &function elevation of T, at some latitude might be a reasonable idealization of the local variation of T, due to enhanced dust/solar absorption in a restricted latitude band, a large scale lati-

tudinal variation of T, , representing the effeet of solar forcing of the rest of the atmosphere, must also be included, if the model is purported to describe an idealized Martian dust storm. The steady symmetric nearly inviscid responses to the background solar forcing and the b-function

319

MARTIAN GREAT DUST STORMS MODELS

-.7

-.6

-.5

-.4

I

-.3

I

- 2

I

-.4

-1

I

-.3

0

I

-.2

I

I1

.2

.3

i

1

4

I,

.5

.6

.7

I

I

I

I

I

.I

.2

.3

.4

.5

-.6

-.5

FIG.

6. Same as Fig. 5, for c = lo-’

forcing are not superposable due to the nonlinearity of the governing equations, but this nonlinearity (in the form of angular momentum conservation) is the essential ingredient that produces the interesting behavior in the pure a-function response problem. The 8 function plus background

-.l

I

.I

set-‘.

solar forcing problem has not been thoroughly analyzed as yet. However, once the characteristics of the solar-forced solutions are described, it will be possible to discuss certain features of the combined problem. In this subsection, the angular momentum conserving response and the response of

EDWIN K. SCHNEIDER

320

0

I

I

I

I

I

I

I

I

I

1

.2

.3

.4

.5

.6

.7

.6

9

I 1.0

Yo -

FIG. 7. Maximum zonal friction velocity, Iu.I,,,, for the response to an isolated heat source of constant strength as a function of the sine (latitude) of the heat source, yO.

the vertically integrated numerical model to T,‘s chosen to represent Martian equinox and solstice conditions will be described. Then it will be possible to describe the combined behavior in certain special cases. Representation of solar forcing. The Leovy and Mintz (1969) calculations showed that the Martian atmosphere is of such small heat capacity that the atmospheric heat fluxes are unimportant for the column-integrated heat balance, but atmospheric heat fluxes are large enough to significantly affect atmospheric temperatures. The zonally averaged net incoming solar radiation and the outgoing thermal radiation at the top of the atmosphere were approximately equal, except over the ice caps. For an optically thin atmosphere in radiative equilibrium, the ground temperature will approach the black body temperature in equilibrium with the net incoming solar flux. This is not far from the case for the Leovy and Mintz (1969) calculations. Accordingly, the radiative equilibrium temperature at the ground is taken to be the ground temperature of the GCM calculations (Leovy and Mintz, 1969, for equinox and southern hemisphere summer and Pollack et al., 1981, for northern hemisphere summer). The solar forcing as represented by Fe in (40) is the vertical mean atmospheric radiative equilibrium temperature, or the atmospheric radiative equilibrium

temperature at approximately H/2. The radiative-convective-condensation equilibrium temperature calculations of Haberle et al. (1982) show that for the clear Martian atmosphere, the winter midlatitude and polar region lapse rate is much closer to isothermal than the summer hemisphere lapse rate. The T, for the solstice simulations is then chosen to be the surface temperature of the GCM calculations extrapolated to Hl 2 using a lapse rate which is -3.O”K km-’ (dust-free atmosphere) for the summer hemisphere and winter hemisphere to 15” latitude, O”K km-’ over the winter ice cap, and varying linearly in y from -3 to 0°K km-* between 15” and the ice cap latitude in the winter hemisphere. The effect of CO2 condensation is taken into account by choosing the minimum value for T, of 150°K. The forms of Fe found for H = 11.6 km are shown in Fig. 8. Angular momentum-conserving response to equinox and solstice Martian solar forcing, The equinox Fe is approximately T, = 133 + 65V’?=-7

(49)

in degrees K. The angular momentum conserving response to this 7, may be found from the small angle approximation expression of Held and Hou (1980). Inserting the relevant numbers, including H = 11.6 km, into their formula, yields Hadley cells extending from

321

MARTIAN GREAT DUST STORMS MODELS meridional cell. That is,

-1 5 y % 1, (50)

u(y,H) = -aam, 24Ok

and from (12) 1 aT -- = Toay

Y/R,

-1 5 y I 1.

(5 1)

Application of the equal area argument, (16), then completes the solution. This solstice-type solution is unusual for a number of reasons. One feature can be illustrated by choosing ,201 -.95

I -.5 SIN

I

I

1

0

5

.95

T &l-+! To

LATITUDE

FIG. 8. Fe chosen to represent solar forcing at equinox (+), southern hemisphere summer(O), and southem hemisphere winter (A). The summer pole in the solstice calculations is taken at y = - 1.

(52)

with pole-to-pole temperature AhTO.Condition (16) then gives

T

-=1-6R+2R

1

Y2

To

the eqUatOr t0 yH = k .45, and WeStedy jets at yn of 49 m set-i. For 1y[ > yn the solution is radiative equilibrium thermal wind balance with no meridional circulation. In order to find the angular momentumconserving solution forced by the solstice T,‘s, note that T, is monotonically decreasing from the summer hemisphere pole to the winter hemisphere pole, and that a i’,/ay is nonzero at the equator. A meridional circulation is necessary whenever the absolute angular momentum at some point of the radiative equilibrium thermal wind is not the absolute angular momentum of the planet at some latitude (Schneider, 1977), and this condition is clearly satisfied when aTJay # 0 at y = 0. Due to the condition that the surface angular momentum is transmitted to the upper streamline at latitude y. by upward motion, T, - T > 0 and m (a T/ay) = 0 at y = y. . Unless y. = 2 1, T has a relative maximum at yo. When aTJay does not change sign (and in some situations where a i’Ja y does change sign), the only angular momentum-conserving solution that satisfies the conditions is y. = ? 1, zero absolute angular momentum zonal winds at 5 = H, and a single pole-to-pole

contrast

(53)

which is inconsistent with the condition T, - 2” > 0 at y = -1 unless R > 2/(3AJ or equivalently H>&=j-~-g

2 (flu)2 (54)

For Martian parameters, H, = (lo/Ah) km, or about 41 km for the southern summer T, and 66 km for the northern summer ?‘, shown in Fig. 8. The condition (54) for the existence of a solstice-type nearly inviscid solution is an analog for symmetric instability in the nearly inviscid model. This condition may be rewritten as (55) where A, is (l/To) times the vertical potential temperature difference from the ground to H, and Ri is the Richardson’s number at the equator, the latitude of largest (au/a& Ri

=

d&W -. VW2 H2

Condition (55) is a stronger or weaker con-

322 straint than

EDWIN K. SCHNEIDER

usual criterion Ri > 1, depending

on whether ($)A, < A, or ($)A,, > A,. If H < H, for some reason, no steady nearly inviscid symmetric circulation exists. A steady symmetric circulation requires friction in this event. In a numerical model, the largest possible His determined by the choice of grid-point heights. For solstice-type forcing, this choice may force H < H, and require a frictionally dominated meridional circulation, where the amount of friction required is related to H. A drag-type friction, as in (41), will alter the solstice-type nearly inviscid solution by adding absolute angular momentum along the upper streamline throughout the summer hemisphere, producing westerlies and a relative maximum of Tin the winter hemisphere, and decreasing the magnitude of ai’/‘/ay.The determination of the correct choice of H for the nearly inviscid model is open. In view of the preceding discussion, a reasonable numerical simulation of the Martian solstice general circulation may require sufficient vertical resolution to represent a very deep troposphere. Response of the vertically integrated numerical model. The numerical model, with

T, chosen to represent the solstice and equinox distributions of incoming solar radiation on Mars, was integrated. The minimum friction criterion was used to determine c. Choices of the parameters, T, (T, and H were made in order to compare the performance of the model to the clear atmosphere GCM calculations of Leovy and Mintz (1969) and Pollack et al. (1981). Another set of calculations with o and H chosen to represent the Martian atmosphere in normal slightly dusty conditions was performed. The choice of H for comparison of the integrated to the GCM calculations is dictated by the GCM vertical grid structure. A choice of H = 1.16 x lo4 m, approximately the height of the highest level of the Leovy and Mintz (1969) calculations at which u is predicted, was used for the comparison cal-

culations. The value for the static stability of v = 1.5 x 10e3 “K m-l (approximately one-third the static stability of an isothermal atmosphere) was taken by inspection of the GCM calculations. The radiative time constant used is 4 Martian days, r = 3.551 x lo5 sec. This value of 7 is given by Haberle et al. (1982) for the clear Martian atmosphere. The model equations (12), (40), (41), and (46) may be nondimensionalized. By the choices of scalings u - Ra, m - H, v (aTO)/( c - T&HuT),and Use- (fIaTo)/ (UT), the resulting nondimensional equations contain only the dimensionless parameters R = gH/(Ra)*, the nondimensional vertical structure constant m, and the nondimensionalized value of the friction coefficient c. Given fixed vertical structure and dimensionless forcing (T,IT”), and using the minimum friction criterion (c = c), the solution to the nondimensional equations depend only on R. Thus, the nearly inviscid numerical model has the same properties as the model of Section 2. The structures and magnitudes of u(y,H) and T do not depend on o and T. The meridional structures of the meridional velocity and surface stresses are also independent of (+ and T, while their magnitudes vary as (VT))‘. The results of the GCM comparison calculations are shown in Fig. 9. For all three of the cases, equinox, southern hemisphere summer, and southern hemisphere winter, the position of the westerly zonal wind maxima (the jet streams) in the GCM calculations is well simulated by the integrated model. Additionally, the relative variation of the speed of the jet stream is similar in the GCM and integrated model calculations, with the magnitudes of the jet stream zonal velocity about 25% larger in the integrated model. Thus, the large differences in the latitude and strength of the Martin GCM jet streams are given by the overturning meridional circulation response to the seasonally varying solar forcing; there is no need to appeal to the nonzero mass flux

MARTIAN

323

GREAT DUST STORMS MODELS

a

b

-lOOl -10

-05 SIN

1201 -0.95

I -05 SIN

0 LATITUDE

0.5

I 0

I 0.5

10 SIN

LATITUDE

SIN

LATITUDE

I 095

LATITUDE

FIG. 9. Results of vertically integrated model response to seasonal variation of solar forcing on Mars for clear atmosphere cases. Equinox results marked with +, southern summer with 0, and southern winter with A. For purposes of comparison, summer pole is at y = - 1. (a) Zonal winds at the top of the atmosphere. (b) Meridional winds at the top of the atmosphere. (c) Atmospheric vertical mean temperature. (d) Surface friction velocity component in the zonal direction.

condensation flow to explain the seasonal variations of the GCM jet streams. Surface stresses in the polar regions will of course be strongly affected by the condensation flow. Large differences in the strength of the tropical and summer hemisphere easterlies between the solstice integrated model and GCM calculations are apparent with the integrated model having much stronger easterlies. Haberle et al. (1982) obtain much

closer agreement between a low vertical resolution two-dimensional model and the Pollack et al. (1981) GCM results in these regions. Apparently, the neglect of eddy momentum flux divergences is then not the dominant cause of the discrepancies. The cause appears to be due to the difference in the form of the zonal momentum equation between the layer representation adopted here and the level representation of the vertically finite differenced models. The inte-

EDWIN K. SCHNEIDER

grated model uses the zonal momentum equation at 5 = H where w = 0. Hence, there is no contribution from vertical advections of zonal momentum. The twolevel-type model applies the zonal momentum equation at an upper level which is midway between the top (w = 0) and the middle of the atmosphere, and the vertical shear of the zonal wind at this level is taken as the vertical shear at the middle level. The vertical advection of zonal momentum by a zonally symmetric circulation acts as a westerly torque on the upper level zonal winds in regions of upward motion and easterly shear with height, and regions of downward motion and westerly shear with height, leading to weaker easterlies and stronger westerlies in the two level than the integrated model. If the vertical advection of zonal momentum is included in the zonal momentum equation (a level representation), results closer to the upper level u of Haberle et al. (1982) are obtained for southern hemisphere winter (equatorial easterlies reduced by a factor of 2). The latitude of the jet stream remains the same with the layer and level representations of the zonal momentum equation, and the speed of the jet stream increases in the level representation, as expected from the previous discussion. As vertical advection of zonal momentum by the mean meridional circulation cannot alter the conservation property of absolute angular momentum, its pseudofrictional effect in low vertical resolution models must be considered to be a liability. The results for u(y,H), VeH’HO(the meridional velocity at the top of the model atmosphere), T, and U*for the three GCM comparison calculations are shown in Figs. 9a-d. The zonal winds (Fig. 9a) show the large response to the seasonal variation of the solar forcing described above. The meridional wind at the top (Fig. 9b) also shows an extremely strong seasonal response to the solar forcing. The solstice meridional winds reach speeds more than an order of magnitude larger than the equinox meridional winds, and the southern summer merid-

ional winds are approximately double those obtained for southern winter. This seasonal variation in the strength of the meridional circulation is qualitatively similar to the results from the GCM calculations, where the maximum upper level meridional winds were approximately 1 m set-’ at equinox, 5 m set-i in southern winter, and 10 m set’ in southern summer. Solar forcing asymmetric about the equator produces a much stronger meridional circulation than forcing symmetric about the equator. The meridional velocities are somewhat stronger at the solstices than those found in the GCM calculations, while the equinox meridional velocities are similar to those produced by the GCM’s. The southern winter meridional wind speed is close to that produced by Haberle et al. (1982) in their GCM comparison calculation. The solstice meridional circulations produced by the integrated model extend to close to the summer hemisphere pole, while in the GCM and two-dimensional models, the summer to winter hemisphere meridional cell begins at 30 to 40” latitude in the summer hemisphere. This difference is due somewhat to the choice of T, in the one-layer model to reflect pure radiative equilibrium rather than radiative-convective equilibrium. Convection in both the GCM and two-dimensional simulations heats the atmosphere most strongly near the subsolar latitude, which would lead to a radiative convective T, with a maximum near the subsolar latitude and a decrease towards the summar hemisphere pole. 7, was chosen to represent the external forcing only in the one-layer calculations, as this is the simplest and most direct choice. The model zonal mean temperatures (Fig. SC) show the features expected from the discussion of Section 3. At equinox, T - T, everywhere, but the small difference between 7 and T, in the region equatorward of the jet streams is sufficient to cause a large difference in this region between the dynamically produced zonal winds and the radiative equilibrium thermal winds (ap-

MARTIAN GREAT DUST STORMS MODELS proximately solid body rotation). At the solstices there is a relative minimum in T at the equator and relative maximum in 7 in the winter hemisphere. The winter hemisphere maximum in T, which is closely connected to the jet stream, becomes more pronounced and moves poleward as the strength of the thermal driving increases, and appears only when the forcing is asymmetric about the equator. The friction velocity, u*, given by (46), is shown for the model calculations in Fig. 9d. The sign of the surface zonal wind is given by the sign of u*. u* is strongly seasonally dependent, with negligible 1u*l at equinox compared to the solstices. Additionally, the largest 1~~1’s for southern summer forcing are almost double those for southern winter forcing. The magnitude of the friction velocity is somewhat larger than that calculated by Haberle ef al. (1982). The zonal surface stresses, -u*l u+I, given by the model agree in relative magnitude for the three seasons with the zonal surface stresses produced by the overturning meridional circulation in the GCM calculations. The results are suggestive as to a mean meridional circulation explanation for the seasonal behavior of Martian great dust storms. The southern summer axially symmetric response to solar forcing is significantly closer to exceeding the threshold friction velocity for saltation, whatever it is, than the southern winter response. The equinox response is much farther from exceeding the threshold friction velocity than either solstice season. The meridional component of the friction velocity is not calculated in the model. This component is expected to be strongest in the tropics at solstice, where the zonal component is weakest, leading to a less extreme latitudinal variation of the total friction velocity than given by the zonal component alone. (u+( maxima at solstice are found in the subtropics of both hemispheres and near the jet stream. A similar structure was found in the diagnostic study of Leovy et al. (1973). Haberle et al. (1982) found a maximum in

325

the surface wind speed of the summer hemisphere at 30” for small and moderate atmospheric dustiness. The subtropics and jet stream regions would then, at first guess, seem to be the most likely regions for dust raising in an axially symmetric atmosphere. Both the tilt of the Martian rotation axis and the distance from sun would appear to have significant effects on 1u* 1, with the obliquity effect being the more important. The angular momentum conserving equinox solution is in good qualitative agreement with the equinox GCM and vertically integrated solutions. At the solstices the minimum friction vertically integrated model and the GCM also are in qualitative agreement, with the primary disagreement resulting from the neglect of w(aula{) pseudofriction in the vertically integrated model. At the solstices, the GCM vertical structure dictates a choice H < H,, so that an angular momentum-conserving solution does not exist by (54). Internal friction has a significant effect on the GCM calculation. Whether this effect is a representation of some real physical process, or an artificial effect necessitated by the choice of the vertical grid structure is a problem that needs to be investigated further. Even the model of Haberle et al. (1982), which extends to H - 35 km is not deep enough to allow minimization of unresolved frictional effects. Insufficient model depth may explain the necessity in their calculations for a critical Richardson’s number criterion Ri, = 2.5 > 1 for numerical stability. Due probably to the presence of suspended dust year round (Pollack et al., 1979), the static stability of the ambient Martian atmosphere is considerably larger than that produced by the clear atmosphere GCM calculations. The height of the Martian troposphere appears to be increased by dust/solar absorption (Haberle et al., 1982). No Martian GCM calculations including the observed atmospheric dustiness have been published. The parameters of the model were adjusted to estimate the change in the axially symmetric response with a more re-

326

EDWIN K. SCHNEIDER

alistic o and H than used in the GCM comparison calculations. As a crude representation of the effects of the small amount of dust, the model static stability was doubled to 3°K km-i and the height of the model was doubled to H = 23.2 km. The same method of extrapolation of the surface temperature to H/2 gave the same Fe’s as used in the previous calculations. The radiative time constant is presumably slightly shortened by dust/solar absorption, but T = 3.351 x lo5 set was taken as before. Calculations were again performed for the minimum internal friction. It is seen from the nondimensional model equations that season by season differences in u(y,H) and F between these results and the previous seasonal calculations are due only to the change (doubling) of the parameter R = gH/(Ra)Z. Changes in the meridional velocity and U* are due to both the increase in R and the static stability (T, while the effect of a change in the radiative time constant is given by the scaling u**, V - 7-l. The results for uo),H) show poleward displacement and strengthening of the jet streams relative to the clear atmosphere results. ,Stronger easterlies are also produced in the solstice seasons (- 120 m/set near the equator for southern hemisphere summer vs -80 m/set). The zonal winds at equinox in upper latitudes are double those of the clear equinox calculation, a consequence of the doubling of the vertical scale, since the zonal winds are in radiative equilibrium thermal wind balance poleward of the jet streams. The meridional winds at the top of the atmosphere are weaker than previously, and at solstice the meridional wind maximum is displaced towards the winter hemisphere pole. The net meridional mass flux, VH, decreases by a factor of 2 in the solstice seasons. This weakening of the meridional circulation is primarily due to the doubling of the static stability and the vertical structure assumptions (ucw~‘~o,o, and 7 = constants with height) used in deriving the

integrated representation of the heat equation. The same meridional mass flux produces a larger adiabatic warming for the same stability when the vertical scale is expanded, due to the decrease of density with height. Also, the meridional circulation is weakened as H + H, as less friction is required for the existence of a steady solution. The temperature maxima in the winter hemisphere dusty atmosphere are stronger and more poleward. The winter hemisphere extratropics are somewhat warmer than in the clear atmosphere case, and the summer hemisphere T is cooler, except for little change near the summer hemisphere pole. The changes as H is increased lead to a closer approach to the angular momentumconserving solstice solution, since a smaller effective friction is needed as H + H,. As H + H,, the winter hemisphere temperature maximum and westerlies approach the winter hemisphere pole, and the westerly jet strengthens. Haberle et al. (1982) possibly could not simulate the winter polar warming observed at the height of the global dust storms due to insufficient model depth. The equinox T is close to that found previously. The zonal component of the friction velocity is decreased by about 20% from the clear atmosphere case, due to the decrease in the meridional mass flux. The relative variation of (u*) for the three seasons is as before, with southern summer having the largest values by a factor of about 2 over southern winter, and a factor of about 10 over the equinox. Hence, speculation on the meridional circulation role in generating global dust storms is as above. Response to 6 function plus solar forcing. In the context of the constant angular

momentum or vertically integrated model, the effect of dust/solar absorption in a narrow latitude band at y. on the steady zonally symmetric atmospheric circulation can be found as the response to a i;, of the form

C?(Y) QY> -=+ NY To To

Yo)

(56)

MARTIAN

GREAT DUST STORMS MODELS

(cf. Eq. (18)). FJy> represents the global scale variation of Fe due to solar forcing of the ambient atmosphere. The analysis of the problem with i’e given by (57) has not been completed. It is possible, however, to make deductions from graphical analysis, as in Section 3. It is convenient to define “local” and “global” responses with reference to the region to which dust added to the atmosphere at latitude y. could be transported. If the meridional cell(s) in communication with the dust at y. could transport it across the equator, the response will be called global. Otherwise the response will be called local. The meridional circulation driven by fe will be considered to have a vertical scale H. The forcing due to the 6 function has a vertical scale Hi, which is the height of the dust column. The distinction as to local or global behavior is made with reference to the implied transport of dust; hence the structure of the meridional circulation at Hi should be the appropriate indicator of “localness” or “globalness.” The simplest case to consider is H1 = H, neglecting modification to the angular momentum-conserving solution by vertical momentum transport at y. when motion along the upper streamline passes from one side of y. to the other. For large enough q, the response is always global no matter what the form of f.,, and this response is essentially identical to the two-celled, u = 0 at y = yo, global response of Section 3. Therefore, if a local response can be shown to exist for a given fe and y. and small enough q, some sort of transition from local to global response must occur as q is increased. If fe is such that the region A affected by the meridional circulation with q = 0, given by -ys 5 y 5 yN, ys, yN > 0, does not reach the poles (typified by the equinox i;,), then a small enough q will always produce a local response when y. is outside of A. Superposition is valid in this case, and there will be a four-celled meridional circulation. If y. is inside of A, the combined response will be a small perturbation to the

327

two-celled q = 0 response for small enough q. The classification of the response will depend on where dust injected at y. would be transported to by the q = 0 solution: global ify, is between the u = 0 rising motion latitude of the q = 0 solution and the equator, and local otherwise. Thus, there is no local response to a summer hemisphere yo for solstice-type Fe,. This result is a cons-quence of the condition of continuity of T. However, dust transport from yo toward the summer hemisphere pole will require the qualitatively different large q-type global response. The situation for which local to global transition is sought in order to claim relevance of this model to the onset of global dust storms is solstice-type solar forcing with y. in the summer hemisphere. In this case, there is no such transition when HI = H. In order to find local behavior for the relevant case, it is necessary to consider the small q solution for H, < H. First, the outer solution with the vertical mean operator applied over H is calculated. Assuming that H is high enough so that internal friction is not necessary, this solution is given above in Section 4 (Eqs. (50)~(51)) (T symmetric about the equator), with q entering into determination of the constants of integration. Next, it is assumed that the b-function forcing induces an angular momentumconserving response in a region ys 5 y 5 yN at HI. The inner solution is thought of as a bubble containing those parcels which pass through the 6 function. The calculation of the structure of this inner solution differs from that of Section 3 in that the exterior (formerly undisturbed) region is not in radiative equilibrium. A local response to the 6 function for solstice-type fe and summer hemisphere y,-, would be bounded meridionally by regions where i;, < Fe (vertical means taken over height HI) in which the vertical velocity is upwards (>O). The vertical velocity in the inner bubble must, of course, be 50 near ys and yN. As T must still be continuous, the inner bubble and exterior region have to be separated by a slop-

328

EDWIN K. SCHNEIDER

ing front at which vertical velocity is discontinuous. Then an upward mean vertical velocity near ys+ and yN- is consistent with the absolute angular momentum at (ys+, Hi) and (YN-, Hi) being the absolute angular momentum at the ground at y. (i.e., the analog of condition (27) does not apply in this case). The end result is that the construction of Section 3 is “superposeable” on the T of the outer solution extrapolated so as to be the mean temperature between the ground and Hr. uCy,Hi) in the outer solution can be approximated by u(y,H) x (HI/~. Redefining q in (57) as (H,/H)q so that the strength of the 6 function integrated over H is independent of H, the construction of Section 3 for the inner solution (before checking for consistency) proceeds as before. Continuity of T at yN and ys gives qualitatively the same result as in Fig. 2: one, two, or three possible solutions with different values of Rq (equivalently H,q). By replacing condition (27) with the condition resulting from the heuristic discussion of Section 3, that the fronts in the consistent solution be statically stable so that the discontinuity in u is relative westerlies to relative easterlies moving away from the equator, the intermediate solution is again rejected. The qualitative behavior of abrupt local to global transition of the dust transporting inner solution is then found when large-scale solstice-type solar forcing is included. The equal area condition (20) with the appropriate Fe(y) replacing TA may be used to determine the strength of the 6 function (H,q), at which a local to global structural transition will occur. Frontal discontinuities, however, may have to be included in the vertical integration of the thermal wind equation for consistency. The effect of this is not known. Since fe > T when q = 0 in the region that will be occupied by the inner solution, it is found (assuming aT/ay = aF//ay> that the width of the region occupied by the inner solution approaches some nonzero value as q + 0 for fixed HI (summer hemisphere yo), and that the width of this limit-

ing region approaches zero as H, + 0. Additionally, there is some value of H, = Hlc < H such that if HI > Hlc , there is no local solution. These results follow from continuity of T, and it is also possible to obtain them from the condition for frontal instability. It therefore seems that a zonally symmetric nearly inviscid model contains the necessary ingredients to simulate an explosive transition between quasisteady local and global dust storm regimes as some critical amount of dust/solar absorption is exceeded. 5. DISCUSSION The steady, zonally symmetric, nearly inviscid model response to thermal forcing has been shown to contain behaviors which may be relevant in understanding the rapid onset of Martian great dust storms. Discontinuous large increases in the scale of the meridional circulation response to enhanced dust/solar absorption in a small latitude band can occur as the amount of absorption is increased a small amount from below to above some critical value. These increases in scale are suggestive of explosive growth of a local to a global dust storm, in terms of dust transport. The model also produces a dynamical response to the seasonal variation of solar insolation that is qualitatively similar to that found in GCM and two-dimensional numerical simulations. The zonally symmetric response (neglecting ice cap mass sources and sinks) has been shown to strongly favor the southern hemisphere summer subtropics for the occurrence of dust raising. The zonally symmetric solar-forced response to timedependent solar forcing would initiate a global dust storm with a local storm in the latitude band and season where global storm precursors are observed to originate for a wide range of possible threshold friction velocities for dust raising. Due to dust/ solar absorption, the local storm would grow slowly in area (and height) until some threshold amount of dust was raised, at

MARTIAN

GREAT DUST STORMS MODELS

which time explosive growth to a global dust storm (by transport) would occur. The global storm would spread dust latitudinally, and significantly alter the atmospheric temperature distribution, with the amount of polar penetration of dust and warming depending on the effective height of the storm. A truly global storm would have to be very deep. The model is, of course, incapable of representing the potentially important effects due to transience (local storm lifetime and diurnal variation of all types) and forcing localized in longitude as well as latitude. In view of the above results, it seems possible that a zonally symmetric numerical atmospheric model including a dust continuity equation could provide a realistic simulation of the present seasonal behavior of Martian great dust storms, without including currently popular tidal or topographic effects. The angular momentum conserving model provides a guide as to how such a model could be designed. The model of Haberle et al. (1982) is a significant first step in the development of such a model. The timing of the occurrence of great dust storms in the proposed model would depend on the choice of the threshold friction velocity for dust raising, 1u*, 1. Presuming that the dependence of 1ulc 1 on surface pressure and temperature is given by some similarity theory, the model would be calibrated by a single constant number. There would be two obvious tests as to the necessity to complicate the model (i.e., with interactive ice caps, tidal motion, topography, or baroclinic instabilities). One test would be the ability of the zonally symmetric model to simulate the observed decay of a great dust storm. The zonally symmetric model must accurately treat the interaction of dust with solar and thermal radiation to perform this test. The other test would be the ability of the model to simulate interannual variability of the occurrence of great dust storms. As the radiative time constants are so short, the meridional circulation would always be in quasiequi-

329

librium with the solar forcing, aside from possible symmetric instability, except during the growth of a great dust storm. Thus, interannual variability in the occurrence of great dust storms would require a mechanism for interannual variability in the radiative forcing of the atmosphere. Two possible ways that a zonally symmetric model could see such variability are very long time scales for dust fallout, and memory of the preceding dust storm(s) contained in dust deposition or extent of the polar ice caps. Finally, the implications of the angular momentum-conserving or minimum friction model for the dependence of the occurrence of great dust storms on Martian orbital parameters will be mentioned. The atmospheric mass (PO) is expected to have varied significantly, perhaps by more than an order of magnitude, as the orbital parameters have changed (Ward et al., 1974). Various theories predict that the threshold friction velocity should vary as pop213(Ward et al., 1974) or approximately as p0-1’2 (Pollack et al., 1976). By assuming that the dynamically induced friction velocity is approximately independent of po, the case has been made (Ward et al., 1974; Pollack, 1979; Pollack and Toon, 1982) that at times of low surface pressure (i.e., low obliquity) global dust storms should disappear, and at times of high surface pressure (i.e., high obliquity) there may be global dust storms year round. The angular momentum-conserving model makes a similar though more detailed prediction, but for entirely different reasons. The angular momentum-conserving model says that for a fixed distribution of solar insolation and H, the zonal winds will be independent of variation of the other parameters (i.e., surface pressure), while the meridional velocities will vary as (UT)-‘. If the static stability, (+, is fixed, (44) then implies that the friction velocity due to the zonally symmetric meridional circulation varies as T-I/~. As the radiative time constant, T, is expected to vary proportional to

330

EDWIN K. SCHNEIDER

po (Gierarch et al., 1970), the dynamically induced friction velocity will vary as po-1’2. Then the ratio of the threshold friction velocity to the dynamically induced friction velocity will vary as po-1’6 or poo depending on which theory for the threshold friction velocity is used, and the variation of p. is a weak effect. Thus, the angular momentumconserving or minimum friction models predict that to a first approximation the variation in po can be neglected in deducing the effect of variations of the orbital parameters. The occurrence of global dust storms then depends only on the instantaneous distribution and intensity of the incoming solar radiation. The constant 7, constant 1ye?1,angular momentum-conserving or minimum friction model calibrated to behave correctly for the present p. and seasonal cycle of solar insolation will also predict the occurrence of great dust storms in past climatic regimes merely by choosing the appropriate T,. By examination of Fig. 9d, it is qualitatively seen that during times of low obliquity (demarked by some critical value) there will be no global dust storms, the same prediction as when dynamical variations are neglected. During periods of high obliquity there may be two seasons (solstices) of global dust storm activity, but global dust storms will never occur in equinox seasons as has been hypothesized. Eccentricity and time of perihelion will also play significant roles when the obliquity is large enough. Quantitative calculations of the year by year occurence and duration of global dust storm seasons in the zonally symmetric model with the orbital parameters changing as they are presumed to have changed would be simple to perform. Comparison of these calculations with the record of layering in the polar caps would be of interest.

ACKNOWLEDGMENT This research has been supported by the National Aeronautics and Space Administration under Grant

NGL-22-007-228. Useful suggestions were made by Dr. R. S. Lindzen and anonymous reviewers.

REFERENCES BINDER, A. B., R. E. ARVIDSON, E. A. GUINNESS, K.

L. JONES, E. C. MORRIS,T. A. MUTCH, D. C. PIERI, C. SAGAN (1977). The geology of the Viking lander I site. .I. Geophys. Res. 82, 4439-4451. CONRATH,B. J. (1981). Planetary-scale wave structure in the Martian atmosphere. Icarus, 48, 2739-2748. DICKINSON, R. E. (1971). Analytic model for zonal winds in the tropics. I. Details of the model and simulation of the gross features of the zonal mean troposphere. Mon. Weather Rev. 99, 501-510. GIERASCH,P. J. (1975). Meridional circulation and the maintenance of the Venus atmospheric rotation. J. Atmos. Sci. 32, 1038-1044. GIERASCH,P. J., AND R. M. GOODY (1973). A model of a Martian Great dust storm. J. Atmos. Sci. 30, 169-179. GIERASCH,P. J., R. M. GOODY, AND P. STONE (1970). The energy balance of planetary atmospheres. Geophys. Fluid Dynam. 1, I-18. HABERLE, R. M., C. B. LEOVY, AND J. M. POLLACK (1982). Some effects of global dust storms on the atmospheric circulation of Mars. Icarus, 50, 322367. HELD, I. M., AND A. Y. Hou (1980). Nonlinear axially symmetric circulations in a nearly inviscid atmosphere. J. Armos. Sci. 37, 515-533. HOU, A. Y. (1981). Thermal Convection in Planetary Atmospheres. Ph.D. thesis, Harvard University, Cambridge, Mass. LEOVY, C. B. (1964). Simple models of thermally driven mesospheric circulations. J. Atmos. Sci. 21, 327-341. LEOVY, C. B. (1979). Martian meteorology. Annu. AND

Rev. Asrron.

Astrophys.

17, 387-413.

LEOVY,C. B., ANDY. MINTZ (1969). Numerical simulation of the atmospheric circulation and climate of Mars. J. Armos. Sci. 26, 1167-I 190. LEOVY, C. B., R. W. ZUREK, AND J. B. POLLACK (1973). Mechanisms for Mars dust storms. J. Atmos. Sci. 30, 749-762.

MARGULES,M. (1906). iiber Temperaturschichtung in station&r bewegter und in ruhender Luft. HannBond. Meteorol.

Z. 243-254.

MARTIN, L. J. (1974). The major Martian dust storms of 1971 and 1973. Icarus 23, 108-115. MARTIN, T. Q., AND H. H. KIEFFER (1979). Thermal infrared properties of the Martian atmosphere 2. The 15 micron band measurements. J. Geophys. Res. 84, 2843-2852. MUTCH, T. A., R. E. ARVIDSON,A. B. BINDER, E. A. GUINNESS, AND E. C. MORRIS(1977). The geology of the Viking lander 2 site. J. Geophys. Res. 82, 4452-4467.

MARTIAN

GREAT DUST STORMS MODELS

POLLACK,J. B. (1979). Climatic change on the terres-

trial planets. Icarus 37, 479-553. POLLACK,J. B., D. S. COLBURN, F. M. FLASER, R. KAHN, C. E. CARLSON,AND D. PIDEK (1979). Properties and effects of dust particles suspended in the Martian atmosphere. J. Geophys. Res. 84, 29292945. POLLACK,J. B., R. HABERLE,R. GREELEY,AND J. IVERSON(1976). Estimates of wind speeds required for particle motion on Mars. Icarus 29, 395-417. POLLACK,J. B., C. B. LEOVY, P. W. GREIMAN,AND Y. MINTZ (1981). A Martian general circulationexperiment with large topography. J. Atmos. Sci. 38,

3-29. POLLACK,J. B., AND 0. B. TOON (1982). Quasi-peri-

331

odic climate change on Mars: A review. Icarus 50, 259-287. SAGAN, C., AND R. A. BAGNOLD(1975). Fluid transport on Earth and aeolian transporton Mars. Icarus 26, 209-218. SCHNEIDER,E. K. (1977). Axially symmetric steadystate models of the basic state for instability and climate studies. Part II. Nonlinear calculations. J.

Atmos. Sci. 34, 280-296. WARD, W. R., B. C. MURRAY, AND M. C. MALIN (1974). Climatic variation on Mars. 2. Evolution of carbon dioxide atmosphere and polar caps. J.

Geophys. Res. 79, 3387-3395. ZUREK, R. W. (1976). Diurnal tide in the Martian at-

mosphere. J. Afmos. Sci. 33, 321-337.