Introduction

In trying to understand ice oscillation and related climatic change, one may first ask: is ice advance/retreat caused by external or internal forcings? The external forcings are agencies outside the Earth system, such as astronomically caused changes in the Earth’s orbital parameters.

Internal forcings, on the other hand, are processes and mechanisms present within the Earth system, giving rise to self-excited glacial fluctuations. The temperature-albedo feedback mechanism associated with long and short wave radiation fluxes at the Earth’s surface are given considerable attention in glacial and climate models (Källenand others, Reference Källén, Crafoord and Ghil1979; Oerlemans and van der Veen, Reference Oerlemans and van der Veen1984). However, the role of the hydrological cycle on the ice advance/retreat is given less attention although it is realized that ice freeze/melt and evaporation/precipitation are contributors to the change of air and ocean temperatures, which in turn influence the ice. The purpose of this paper is to investigate, in its simplest form, a possible feedback mechanism in the air—ice-ocean system regarding the hydrological cycle. The ice advance/ retreat on the continent is not included in this work.

Figure 1 shows the main physical processes in the coupled system. From the global point of view, ice advance isolates the ocean from the atmosphere and reduces the heat loss from the ocean surface, which in turn decreases the air temperature and increases the ocean temperature. The air and ocean feed back into the ice through two different mechanisms (Fig. 2): (1) thermodynamical feedback; the increase of ocean temperature melts the ice and causes ice retreat (negative feedback between ice and ocean), however, the decrease of air temperature causes freezing and makes the ice further advance (positive feedback between ice and air): (2) mechanical feedback; the modification of the atmospheric and oceanic temperatures varies the surface evaporation rate, which in turn changes the sea-level height and the ice accumulation rate. Both effects lead to a further change of the ice flow.

Fig. 1. Physical processes in the coupled air—ice-ocean system.

Model Description

Ice advance/retreat

The total surface, ice covered, and open ocean areas are set to be N _i, and 1 − N_i respectively (Fig. 1). The ice is assumed to be isotropic with thickness H_i Following Weertman (Reference Weertman1957), the ice-spreading rate is computed by:

where ρ_i is the ice density. H₁ is the height between the ice top and the sea-level (Fig. 1). Values of n vary from about 1.5 to 4.2 with a mean of about 3, and for randomly oriented polycrystalline ice at −10°C and n = 3, a value A = 3 × 10⁻⁸a⁻¹ k Pa⁻³ is reasonable for ice stress (Weertman, Reference Weertman1973). The angle brackets denote values averaged over ice thickness. The time rate of change of ice coverage N_i contains mechanical and thermodynamical (ice edge freezing/melting) parts:

Fig. 2. Feedback mechanisms among air, ice, and ocean.

where H_i is the ice thickness, and F is the freezing rate at the ice edge. Following Killworth (Reference Killworth1979), the heat flux is accomplished through an exchange coefficient Κ and is of the form K(T − T_f), where T_f represents the freezing point of sea water. Therefore the freezing rate is computed by:

where K_wi, K_ai are the heat exchange coefficients for water-ice and air-ice, respectively. L_i is the latent heat of ice. If T_w, T_a are greater than T_f , F is taken as a negative value, which means ice melting.

Time rate of change of ice thickness

Based on the equation of continuity, Shumskiy (Reference Shumskiy1965) proposed a method of determining ice-thickness changes by comparing the accumulation rate with horizontal strain-rates. The time rate of change of ice thickness becomes:

where c_s is ice accumulation rate, and V_i is vertical mean ice velocity, which links to ice thickness and bedrock topography. It is commonly believed that the non-linear inertial term –V _i -∇H_i in ice-sheet models relates to the long time scale of glacial cycles; however, for simplicity this term is neglected from the present model, i.e.:

which means that the current coupled system excludes the variations with long time scales comparable to the glacial cycles.

Time rate of change of ocean depth

From a thermodynamical point of view, the surface evaporation and ice freezing/melting are two major processes to change the ocean elevation, i.e.:

where H_w is the depth of the ocean, ρ_w is the sea-water density, and E is the surface evaporation rate computed by the bulk formula:

where C _D is the drag coefficient, V _a is the surface wind, q_s is the saturated mixing ratio, and q_a is mixing ratio near the ocean surface.

Heat equations for air and ocean

The heat equations for the atmosphere and ocean are:

where T_a, T_w are temperatures of air and ocean. H_a is the depth of the atmosphere. Q_ai , Q_wi are, respectively, the heat fluxes from air and ocean to ice and Q_wa is the heat flux from ocean to atmosphere. Equations (2), (4), (5), (6), and (7) are the basic equations of our model.

Mean State and Perturbations

Parameterization of perturbations

Mean state of the coupled system (N_i, H_i,H_w, T_a, T_w, E) is obtained by the zero time rate of change in the basic Equations (2), (4), (5), (6) and (7). When the coupled system is perturbed from its mean state, the feedback mechanisms make the perturbation either grow or dampen. Assume that the perturbation of the ocean-surface evaporation rate is caused by the fluctuations of sea-surface temperature (T′_w ), surface wind speed (V′ _a), and the ice coverage (N′ _i ), i.e.:

Ice advance results in an increase of the ocean temperature and a decrease of the surface air temperature. The greater ice coverage will partly insulate the atmosphere from the ocean, decreasing the evaporation (see term (3)); but in locations where the ocean and atmosphere are in direct contact, the warmer ocean temperatures will increase the ocean evaporation rate, which is called the local SST-evaporation coupling (see term (1)). Also the colder air temperatures will strengthen the thermally driven surface wind (due to the larger surface air temperature gradient) and consequently increase the evaporation rate, which is called the wind-evaporation coupling (see term (2)). The net effect as to whether the ice advance increases or decreases the evaporation rate is thus uncertain. According to Chu’s (Reference Chu1989) thermally driven surface-wind formula, the term (2) in Equation (8a) is an order of magnitude smaller than the term (1), and therefore is neglected from the equation, leading to:

where γ(=C _D|V _a|) is the moisture exchange coefficient.

The relationships among the surface evaporation rate E′, the precipitation rate on the ocean

, and the ice accumulation rate are simplified as:

which means the surface evaporation is the source for precipitation and ice accumulation. Following Killworth (Reference Killworth1979), the heat fluxes are computed by:

where c_pa, c_pw are the specific heats for air and ocean, respectively, and K_ai K_wa , and K_wi are heat-exchange coefficients between the corresponding two components.

The perturbation of lateral freezing rate is computed by:

where L_i is the latent heat of the ice.

Timescales of perturbations

Four main timescales appear in this coupled system. They are:

which is the atmospheric relaxation timescale:

which is the oceanic relaxation timescale,

which is the timescale for the change of evaporation rate due to the ocean warming/cooling induced by the ice advance/retreat, and:

which is the ice-spreading timescale.

By using the values listed in Table I, these four timescales are estimated as:

where we assume that |Q_wa|^>>|Q_wi|.

Since the atmospheric and oceanic relaxation timescales τ _a and τ _w are so much shorter than the other two timescales τ _s and τ _e, the atmospheric and oceanic temperature perturbations

and almost “instant¬aneously” follow the ice advance/retreat processes, i.e.:

where

. Equation (17) clearly shows that the ice advance causes the cooling of the atmosphere and the warming of the ocean . Substitution of Equation (17) into Equation (11) leads to the connection of lateral freezing rate F′ to :

where

.

Stability parameters σ and v

Linearization of the basic Equations (2), (4), and (5) leads to:

where

, and σ is a non-dimensional parameter defined by:

For σ > 0, the ice has a “restoring force” induced by the ice advance/retreat. The mechanical feedback of the air and ocean into the ice can be easily noticed from the feedback of

and change due to , into the ice advance/retreat through the first term in the righthand side of Equation (19). This feedback largely depends on the parameter v ≡ H _w/H ₁. The larger the value of v, the stronger this feedback mechanism will be. Figure 3 shows the relationships among the model variables.

Fig. 3. Relationships among the model variables.

The general solutions of Equations (19), (20) and (21) have the following forms:

where b_j and d_j (j = 1,2,3) are the integral constants, and ω ₁, ω ₂, and ω ₃ are the eigenvalues which are the roots of the third-order algebraic equation:

Results

We compute all roots of Equation (24) for different values of the parameters σ and v, and obtain three roots at each pair of the parameters (σ, v). One root among the three has a negative real part throughout the whole parameter space, representing the damping modes, in which we are not interested here. The other two roots, ω₁ and ω₂, have positive real parts somewhere in the parameter space, representing the existence of growing modes in certain parameter ranges. Figure 4 shows the dependence of growth rate (10⁻¹²s⁻¹) on σ for three different values of v. The horizontal axis σ changes from 0 to 10, representing the increase of thermodynamic stability. Oscillation takes place only in certain range of σ (Table II), i.e. for ω _i ≠ 0. Out of this range, i.e. for σ < σ_L, or σ > σ_H, ω _i = 0, which represents non-oscillatory modes. Figure 5 shows the dependence of the period (2π/ω _i) on σ for three different values of v.

Fig. 4. Dependence of growth rate ω_r on stability parameters σ and v.

Fig. 5. Dependence of period of oscillatory modes, 2π/|ω _i|, on stability parameters σ and v.

The results are listed below:

(1) The two parameters σ and v play different roles in the coupled system. The negative feedback measured by σ stabilizes the system; however, the positive feedback measured by v destabilizes the system. Because of the existence of positive feedback mechanisms, the coupled system is unstable unless the parameter σ exceeds some critical value σ_c. Figure 4 indicates that σ_c = 4.25.

(2) Four different types of modes are generated due to the relative strength of positive to negative mechanism:

(3) The parameter v serves as an instability parameter. In the growing regime (Fig. 4), the larger the value of v, the larger the growth rate will be. The mean growth rate of the growing modes varies from 2–6 × 10⁻¹² s⁻¹. The doubling time during which the perturbation doubles its strength is:

(4) The parameter v also largely impacts on the periodicity of the coupled system. First, v affects the critical values of σ_L and σ_Η. The larger the value of v, the smaller the range of oscillatory mode will be. The mean period varies from 2000–3000 years (for v = 200) to 4000–12 000 years (for v = 100).

Conclusion

This conceptual air-ice-ocean coupled model shows a possible positive/negative feed-back mechanism, induced by the hydrological cycle, among the atmosphere, ice, and ocean. The theory predicts the generation of oscillations in ice coverage, ice thickness, air and ocean temperatures, and sea-level height, on the time scale of 10³–10⁴ years.

Acknowledgements

This research was sponsored by the U.S. National Science Foundation, the U.S. Office of Naval Research, and the Naval Postgraduate School.