## Abstract

Classic climatic models use constitutive laws without any response time. A more realistic approach to the natural processes governing climate dynamics must introduce response time for heat and radiation fluxes. Extended irreversible thermodynamics (EIT) is a good thermodynamical framework for introducing nonclassical constitutive laws. In the present study EIT has been used to analyze a Budyko–Sellers one-dimensional energy-balance model developed by G. R. North. The results present self-sustained periodic oscillations when the response time is greater than a critical value. The high-frequency (few kiloyears) damped and nondamped oscillations obtained can be related to abrupt climatic changes without any variation in the external forcing of the system.

## 1. Introduction

Energy balance models (EBM) developed by Budyko (1969) and Sellers (1969) have been a fruitful field of study for diagnosing climatic behavior due to their simple formulation and the qualitative efficiency of their results. The initial studies were mainly related to sensitivity studies of the steady state of the predicted climate (see, e.g., Schneider and Gal-Chen 1973; Held and Suarez 1974; Chýlek and Coakley 1975). Many authors have also analyzed the time evolution of thermodynamic variables using these types of models including seasonal effects (e.g., North and Coakley 1979; Warren and Schneider 1979). However, the time variation studies of EBM have been done by varying the external forcing of the system. This is because internal forcing is not enough to produce periodic or aperiodic damped and nondamped oscillations related to natural climatic variability. With the aim of solving this limitation, statistical dynamics models were developed (see, e.g., Hasselmann 1976; Lemke 1977; Robock 1978), which take into account a stochastic forcing that simulates the intrinsic variability of the system. Nevertheless, using a deterministic EBM with delay time in the ice–albedo feedback, Bhattacharya et al. (1982) obtained self-sustained oscillations that were quasiperiodic or aperiodic. Recently, thermodynamic climate models have been used for studying the role of the rate of entropy production in the climatic system (e.g., Nicolis and Nicolis 1980; Essex 1984; Li and Chýlek 1994; Li et al. 1994).

These models use classic thermodynamic laws, which represent an immediate adjustment between the fluxes and the forces, that is, implicitly accepting that the thermal signal propagates at infinite velocity. The assumption of delay times for the fluxes would allow a more realistic finite speed of propagation. In fact, infinite speed is not a problem since it is generally recognized that the thermal diffusion approximation, Fourier’s law, is only expected to hold in an average ensemble sense. If the velocity fluctuations have correlation times of the order of a few days, it might take the order of 10 times this to form an ensemble in the ergodic sense. Hence, any result of an atmospheric EBM calculated with a response time lower than a month should be considered suspect. Moreover, even though from a physical point of view and also in order to have a well-based model, it is interesting to explore the consequences of the introduction of delay times in an EBM model. Following the results of Hasselmann (1979) high values of the delay time can be expected if we take into account the characteristic timescale of the deep ocean layer. Oceanic circulation is assumed to play an important role in some climatic events, for example, during the Younger Dryas, when a collapse of the Atlantic thermohaline circulation, which initiates the cold event with an important reduction of deep-water formation, is accepted (Stocker and Wright 1996).

Our study goes beyond the mere introduction of delay times because this assumption implies that the fluxes are considered independent variables. We work outside the classic irreversible thermodynamics (CIT) theory, because the CIT uses constitutive laws without any response time and the hypothesis of local equilibrium. The extended irreversible thermodynamics (EIT) theory postulates a generalized entropy that depends on dissipative fluxes and on classic variables (Jou et al. 1996) and it has been used in many condensed matter systems (e.g., Criado-Sancho and Llebot 1993; Fort and Llebot 1996) with very satisfactory results in the high-frequency range (short timescales) where the CIT cannot be applied.

The purpose of this paper is to show how the introduction of delay times for the fluxes in the context of the EIT theory leads to self-sustained periodic oscillations, which could explain internal high-frequency variability (few kiloyears) present in some climatic events. Hence, we have the possibility of studying other thermodynamic variables in the climatic system (e.g., the rate of entropy production) using a simple one-dimensional model with a larger range of application than the classic EBM.

In this work, we use the Budyko–Sellers one-dimensional model developed by North (1975a,b) (section 2a) where a delay time for the heat flux is introduced (section 2b). An analytical solution for time evolution is obtained by assuming an ice-free earth (section 3a). In this case, we show that the amplitude of the temperature (difference between the pole and equator values) can behave nonmonotonically for certain delay times. We study situations closer to the present-day conditions with ice caps in section 3b, where the solutions are obtained by numerical integration. The results give a periodic behavior of planetary temperature when the response time is higher than a critical value. This critical value is found to depend on the planet’s emissivity and conductivity. An analogy with recent studies of the Younger Dryas and, in general, of abrupt high-frequency events due to internal climatic variability (Mikolajewicz and Maier-Reimer 1990; Paillard and Labeyrie 1994; Stocker and Wright 1996) is introduced in section 4 by using a time-dependent response time. We conclude the paper in section 5 with a discussion of the results obtained and, also, a comparison of these with the solutions obtained by Bhattacharya et al. (1982).

## 2. Model

### a. One-dimensional model of North

The model used is a Budyko–Sellers one-dimensional type. The basic equation of the one-dimensional horizontal nonstochastic model is the first law of thermodynamics:

where *C* is the heat capacity, **R** and **J**_{q} are the radiation and heat fluxes respectively, and *T* is the temperature. These variables are a function of time *t* and of latitude *ϕ*, where *x* = sin*ϕ*.

Following the CIT theory, Fourier’s law reads

where *λ*_{q} is the thermal conductivity of the system. In Eq. (2) the heat of fusion is not included although its effect can be qualitatively understood in the model as a natural tendency to reduce the observed variations of temperature and ice sheets. If one finally decides to take into account the heat of fusion, a rough analysis of its importance on the model leads to a different critical response time whose value increases ∼90 years. To simplify our development, however, latent heat has been omitted, this becoming a constant limitation of the model. A detailed study using a one-dimensional discrete model would solve this problem.

The dimensionless gradient operator is defined as (e.g., Li and Chýlek 1994)

where **a**_{x} is the unit vector toward the pole.

Budyko–Sellers models parameterize the divergence of radiation flux as (Budyko 1969)

where the subscripts *s* and *l* indicate shortwave and longwave radiation, **R** = **R**_{l} + **R**_{s}, and *T* is temperature in degrees Celsius. Linear dependence on the surface temperature for the longwave radiation flux in (5) can be accepted due to the range of temperatures considered and to the effect of atmospheric longwave absorption. Here *A* and *B* are empirical constants. Shortwave radiation *Q*(*x*) is defined as *Q*(*x*) = *QS*(*x*) where *Q* is the solar constant divided by four and *S*(*x*) the latitudinal distribution of the incoming radiation. For the coalbedo *a*(*x*; *T*) (*a* = l − *α* with *α* the albedo), we take the step-function form suggested by North (1975b)

where *x*_{s} is the sine of the latitude of the edge of the ice sheet, defined as the latitude where the temperature is 10°C below zero,

The solution of this model is obtained by following the development carried out by North (1975a,b) expanding the variables of the system in series of Legendre polynomials. We restrict ourselves to the second order, eliminating the first-order polynomial because of the assumption of symmetric hemispheres in the climatic model. This yields

where *S*_{2} and *a*_{2} have negative values due to the behavior of the second-order Legendre polynomial *P*_{2}(*x*) and *T*_{2} represents the latitudinal distribution of temperatures. For a real case it would also be negative. The values for the planetary temperature *T*_{0}, the incoming radiation *Q,* and the coalbedo for an ice-free earth *a*_{0} correspond to a meridional integration of (9)–(11) respectively.

Substituting (9)–(11) in (1) and using the parameterizations for the radiation field (5)–(6), the equations that determine the evolution of planetary temperature and the variation of pole–equator temperature (1.5 times *T*_{2}) are

where *H*_{m}(*x*_{s}) for *m* = 0, 2 is defined as

### b. Response time

The delay time *τ*_{q} for the heat flux is introduced assuming a Maxwell–Cattaneo-type equation

which has been applied in the description of hyperbolic heat transport (Criado-Sancho and Llebot 1993). Furthermore, the assumption of (15) instead of Fourier’s law has given satisfactory results explaining the behavior of systems in the high-frequency range in other areas (e.g., hydrodynamics, rheology, . . . ; Jou et al. 1996). In fact, for any physical theory, it would be desirable to have a finite speed of propagation of thermal signals and Eq. (15) is one of the most frequently used alternatives. However, the finite speed would be also obtained including in (15) second, third, and so on, time derivatives of the dissipative fluxes. Nevertheless, it is shown that the effect of the second and higher derivatives in (15) can be included in a simple Maxwell–Cattaneo equation by using an effective response time (Dedeurwaerdere et al. 1996; Jou et al. 1996).

Also Eq. (15) can be deduced if we assume that the heat flux at time *t* + *τ*_{q} is related to the temperature gradient at time *t,* expanding until the first order,

Note that in the steady state, Eq. (15) is Fourier’s law.

The value for the microscopic delay time depends on the system and takes a unique and fixed value. For example, the relaxation time for a monoatomical gas is (Jou et al. 1996)

where *l* is the mean free path, *λ*^{′}_{q} the conductivity, and *C*_{c} the heat capacity.

For climatic models we assume a form like (17) for the delay time. Then, denoting *k*_{qh} the horizontal thermal diffusivity of the system (*k*_{qh} = *λ*^{′}_{q}/*C*_{c}) and defining the characteristic length *l* in function of the delay time as *l* = *τ*_{q}*U* with *U* the mean meridional velocity of the system analyzed (here, it corresponds to the atmospheric, mixed ocean layer or deep ocean layer velocities), *τ*_{q} takes the form

For the atmospheric circulation *k*_{qh} ∼ 2 × 10^{3} m^{2} s^{−1} with *U* ∼ 1 m s^{−1}, producing a value of *τ*_{qa} ∼ 0.7 h. In reference to oceanic circulation, using *k*_{qh} ∼ 10^{4} m^{2} s^{−1} with *U* ∼ 2.5 cm s^{−1} and *U* ∼ 0.5 mm s^{−1} for the mixed ocean layer and the deep ocean layer, respectively, Eq. (18) gives values of *τ*_{q} in agreement with that proposed by Hasselmann (1979) (*τ*_{qm} ∼ 6 months, *τ*_{qd} ∼ 1000 years). Accepting that the atmospheric behavior is strongly linked to the oceanic circulation, we can assume a response time for the present climatic state as an intermediate value of *τ*_{qa} and *τ*_{qm} because the model analyzed does not differentiate between the two main phenomena of heat transport. Nevertheless, higher values can be admitted if a strong reduction in oceanic circulation is considered (even higher than the current *τ*_{qd}).

In the climate simulation using discrete models, the maximum size (*l*_{max}) of the grid should correspond to the characteristic length in function of the response time used,^{1}

Equation (18) with *τ*_{q} = 0 leads to the classical result of an infinite speed of propagation of thermal disturbances. Equation (18) also shows how the response time can change its value if, for example, the mean meridional velocity *U* varies.

It would be also possible to consider a response time for the radiation flux but in this work we accept a faster adjustment for radiation than for convection, so no response time for the radiation flux is taken into account.

Using (15) instead of Fourier’s law, the equation that governs the evolution of *T*_{2} takes the form

The equation for the evolution of the planetary temperature *T*_{0}, remains the same, Eq. (12), because in North’s EBM diffusive model the heat flux is only a function of *T*_{2} and, in consequence, the effect of the delay time is only incorporated in the equation for its time evolution. Nevertheless, *T*_{0} depends on the response time as (12) and (20) are coupled through the expression for the ice-sheet value, Eq. (8).

In the next section we analyze the solutions obtained for an ice-free earth and for a more realistic case of an earth with ice caps.

## 3. Solutions

### a. Ice-free earth

Equations (12) and (20) can be solved analytically assuming an ice-free earth. In this case, without an ice sheet for any time (*x*_{s} = 1), (12) and (20) take the form

The solution for the planetary temperature is a damping exponential

where *T*_{00} = *T*_{0}(*t* = 0) and *T*_{0∞} = *T*_{0}(*t* = ∞) = (*Qa*_{0} − *A*)/*B.* Thus, the final value in the steady state only depends on the numerical values assumed for the parameters of the system. Here *T*_{0} will always evolve monotonically and it will be convex or concave depending on the value of the initial condition *T*_{00}.

#### 1) Classic solution

The assumption of Fourier’s classic law (*τ*_{q} = 0) instead of Maxwell and Cattaneo’s equation for the heat flux, Eq. (15), gives another monotonical behavior for the amplitude of the temperature *T*_{2}:

where *T*_{20} = *T*_{2}(*t* = 0) and *T*_{2∞} = *T*_{2}(*t* = ∞) = (*Qa*_{0}*S*_{2})/(*B* + 6*λ*_{q}).

The value of the planetary temperature *T*_{0} is found assuming the earth as a whole and, therefore, only taking into account the emitted and incoming radiation, whereas the value of *T*_{2} is obtained taking into account the heat flux. In that case we get a quicker damping for *T*_{2} than for *T*_{0} as we can see by comparing (24) with (23).

#### 2) Response time

The solution of *T*_{2} including the response time takes the form of a double damping exponential

where the dimensionless time *t*′ is defined as *t*′ ≡ *t*/*τ*_{q} and *T*_{q} is the time dependence part of the heat flux in degrees Celsius (**J**_{q}(*t, x*) = −*λ*_{q}*T*_{q}(*t*)**∇***P*_{2}(*x*)). Of course, *T*_{q} has to satisfy the Maxwell–Cattaneo-type equation

and the constants *D, E, D*′, and *E*′ satisfy

*d* and *e* are known and are a function of the numerical values used for the parameters of the system

where *d* < *e,* and

Thus, (27a)–(27c) is a set of three equations and four unknowns that allow us to choose an arbitrary value for the initial heat flux of the system. In the classic solution this possibility does not exist because the heat flux is linked to the temperature gradient by the phenomenological law, Eq. (2). Fixing the classic parameters used in the system, the variations in response time are the only ones that can affect the values of the exponential terms *d* and *e.* From (29a)–(29b) the range of the response time reaching imaginary values for *β* is

In this case, a damped simple harmonic oscillator (SHO) for *T*_{2} and *T*_{q} will be obtained.

With the values of the parameters used in Nicolis and Nicolis (1980), the solution is a damped SHO when

For other values, the system evolves as a damping exponential type solution.

### b. Earth with ice caps

Taking into account the role of the ice caps in the climatic system, the variables *T*_{0}, *T*_{2}, *T*_{q}, and, in consequence, *x*_{s} are obtained by numerical integration using a standard Runge–Kutta method. In accordance with Hasselmann (1979) and Bhattacharya et al. (1982) for an evolution in a long period (10^{3}–10^{5} yr), we assume a value for the heat capacity of the system taking the deep ocean layer into account (*C* ≈ *C*_{d} ∼ 10^{10} J m^{−2} K^{−1}). For a short period (1–10 yr) *C* takes the mixed ocean layer into account (*C*_{m} ∼ 3 × 10^{8} J m^{−2} K^{−1}).

#### 1) Classic solution

With *τ*_{q} = 0, the planetary temperature *T*_{0} and its amplitude *T*_{2} can evolve nonmonotonically, reaching two possible stable steady states: one related to the present state (*T*_{0} = 14.9°C and *T*_{2} = −28.2°C) and a second one to a completely ice-covered earth (*T*_{0} = −53.9°C and *T*_{2} = −12.3°C).

#### 2) Response time

The introduction of the delay time for the heat flux can completely modify the behavior of the system. A self-sustained periodic solution is obtained with a *τ*_{q} higher than a critical value. This solution implies a change in the behavior of the climatic model in the sense of Lorenz’ classification (Peixoto and Oort 1992). Thus, the system can be almost-intransitive or transitive depending on the value used for the response time. The classic solution produces only the transitive behavior.

The nonlinearity produced by linking *T*_{0} and *T*_{2} through Eq. (8), indicating the ice–albedo feedback, is important enough to maintain the oscillations when a certain delay time is accepted. This can be intuitively understood from the time evolution of *T*_{0}, *T*_{2}, and *T*_{q}, Eqs. (12), (13), and (26). If we assume a high value of *τ*_{q} with the condition *τ*_{q}*dT*_{q}/*dt* ≫ *T*_{q}, Eq. (13) takes the form

which is a nonhomogeneous damped SHO equation. For an ice-variable earth, the term on the rhs of (32) is not a constant, and it could balance the contribution of the damping term. In this case, *T*_{q} would be a nondamped SHO. Then, from (26), *T*_{2} would also evolve in the same form. Also *x*_{s} and *T*_{0} would have similar periodic behavior due to the condition of the ice-sheet edge Eq. (8). For an ice-free earth the results show damped oscillations because the constant value in the rhs of (32) does not balance the damping term. Therefore, in order to obtain periodic behavior in a diffusive EBM model the assumption of a response time for the heat flux is not enough. To produce the nondamped SHO in *T*_{0}, it is necessary to include the ice–albedo feedback via Eq. (8), which links the planetary temperature with its amplitude *T*_{2}, directly related to the behavior of *T*_{q}.

In Fig. 1, we can see the time evolution for *T*_{0} and *x*_{s} with a delay time of *τ*_{q} ≈ 300 yr and initial conditions *T*_{00} = 14.9°C, *T*_{20} = *T*_{q0} = −30°C. The system reaches the classical steady state but shows damping oscillations each with a period of approximately 2000 yr.

For higher values of the response time the self-sustained solution is obtained, as we can observe in Fig. 2, by using *τ*_{q} ≈ 500 yr. The amplitude of the oscillations in *T*_{0} is about 2.7°C with a period of 2000 yr. With the same period there is an ice-sheet variation between 64° and 90° lat with roughly 500 year of ice-free earth conditions in each oscillation. In this case, the temperature *T*_{2} and the heat flux also presents oscillating behavior. The values of *T*_{2} are between −23.9°C and −32.5°C. Its period is identical to that of *T*_{0} and *x*_{s}. The evolution of *T*_{2} is delayed in reference to *T*_{0} with a numerical value for this time lag of about response time *τ*_{q}.

These results have a similar period to that presented by Bhattacharya et al. (1982) using a delay time for the albedo of 1000 yr. However, the amplitude is in comparison, higher in our case [2.7°C vs 0.8°C of Bhattacharya et al. (1982)]. Nevertheless, we can obtain the same period in the evolution, but with a lower amplitude, if we choose the initial conditions closer to the values of the steady state obtained in the classical case. Furthermore, for the initial conditions that imply ice-free earth events the same periodic solution is obtained. This can be seen in Fig. 3, which shows the phase-space trajectories of the system for different initial conditions where *T*_{0} and *T*_{2} are the coordinate axes. On the other hand, a higher value of the response time produces an increase in the period and in the amplitude of the oscillations.

The critical value for the response time has been evaluated for different cases by changing the values of global system emissivity *B,* and conductivity *λ*_{q}. The results are shown in Table 1. An increase in the planet’s emissivity leads to a decrease in the surface temperature, producing an increase in the critical value of the response time *τ*_{qc}. On the other hand, an increase in *λ*_{q} produces a decrease in the critical value of the delay time. In fact, the dependence of *τ*_{qc} on *B* and *λ*_{q} can be understood from Eq. (32). Accepting the pulse frequency of the oscillations *ω* to be of the order of *τ*^{−1}_{qc}, we find a dependence of *τ*_{qc} with the parameters of the system as *τ*_{qc} ∼ *C*/(6*λ*_{q}), in accordance with the behavior shown in Table 1, where an increase in conductivity leads to a decrease in *τ*_{qc}. On the other hand, an increase in global planetary emissivity, *B,* produces a higher damping term in (32). This increases the domain of the delay times that leads to a damped SHO. Therefore, a higher *τ*_{qc} is necessary in order to obtain the nondamped SHO solution.

In Table 1 the value of the ice-sheet edge at the steady state (*x*_{se}) is also shown and it is obtained by using the classic assumption (*τ*_{q} = 0). It can be seen as for *x*_{se} closer to ice-free earth conditions (*x*_{s} = 1), the critical value for the response time *τ*_{qc} is lower. Nevertheless, the amplitude of the periodic behavior is smaller when *x*_{se} increases, becoming zero when *x*_{se} = 1.

Moreover, when the response time is much higher than the critical value, a maximum in the planetary temperature, denoted by *T*_{0max}, is obtained. This value depends on the numerical values used for the parameterizations. For example, *T*_{0max}, = 16.07°C when the numerical values are those used in Fig. 1, but with *τ*_{q} > 2200 yr, *T*_{0max} = 18.07°C if we use the planetary emissivity of *B* = 1.4 W m^{−2} °C^{−1} and, finally, *T*_{0max} = 13.83°C for a solar constant of *Q* = 335 W m^{−2}.

## 4. Time-dependent response time

Recent results obtained through box model (Paillard and Labeyrie 1994), and two-dimensional (Stocker and Wright 1996) and three-dimensional (Mikolajewicz and Maier-Reimer 1990) models have been focused on the importance of the oceanic dynamics for producing abrupt changes in short period events (e.g., the Younger Dryas or Heinrich events). In these cases, the value of the delay time assumed in our work can be related to the variation of the thermohaline circulation in the North Atlantic Ocean according with the assumption of an important addition of freshwater. Therefore, in a real case, the response time would be time dependent. Then, we can accept a Gaussian form for the delay time as

where *σ* indicates the half-width of the distribution of the delay time assuming a maximum value for *τ*_{q} when the oceanic thermohaline circulation was lower, following the idea of Eq. (18).

We show in Fig. 4, the results obtained using a Gaussian form for the response time like (33) with *τ*_{q0} = 2500 yr (deep ocean layer) and *σ* = 900 yr centered at 4000 yr after the beginning, as can be observed in the graph. The minimum value of the delay time has been assumed to be 6 months in accordance with the values for the mixed ocean layer circulation. We can see the variation in almost 4°C in the global temperature over 3000 yr. The initial variation of the ice-sheet edge is found to be between 74° and 58° latitude with a maximum value of 77° before reaching the steady state (73.7°). This behavior would satisfy some abrupt climatic change events of high frequency (e.g., the Younger Dryas). However, it can be also related to events of low variation in the temperature because a decrease in *σ* leads to a decrease in the amplitude of the peak shown in Fig 4. Low amplitude in the oscillations is also expected if the initial conditions of the system are not far from the stationary state of the classic case (*τ*_{q} = 0), or if *τ*_{q0} takes a lower value.

## 5. Conclusions

The internal variability in the climate system is introduced in a simple EBM North model through the assumption of a delay time for the heat flux. The model, therefore, has two independent thermodynamic variables, the temperature and the heat flux with their corresponding evolution equations: the first law of thermodynamics and a Maxwell–Cattaneo-type equation, respectively. The Maxwell–Cattaneo equation for the heat flux is the extension of Fourier’s law into the high-frequency range and it has been used in the framework of extended irreversible thermodynamics (EIT) in order to obtain a finite speed of propagation of thermal disturbances. Other phenomenological laws could overcome this problem but, as Jou et al. (1996) have shown, it would produce slight modifications in the results obtained assuming the Maxwell–Cattaneo equation. The delay time has been calculated through the thermal diffusivity of the system and with an average value for the mean meridional velocity of the subsystem being considered. Hence, values for the atmospheric, mixed ocean, and deep ocean response times have been evaluated.

Although this study does not consider the dynamics of the system, the assumption of delay time introduces dynamical considerations in the simplest way, without going out of the context of a thermodynamics theory (EIT). This fact allows us to maintain the possibility of studying other thermodynamic variables of climatic interest (e.g., entropy production) by using the powerful techniques developed for this type of theory. Moreover, we have pointed out the existence of a numerical response time related to the grid spacing of the model. Thus, the introduction of delay times in a general climate model should not only take the physical description into account but also a computational contribution, this being more important when the grid spacing of the model becomes rougher. In the present paper, we have not considered the computational delay time because in North’s model the variables are continuous in latitude following the behavior of Legendre polynomials.

We have compared the classic solution versus the solution obtained with response time, analyzing two different cases. For an ice-free earth, there is a decoupling between the two components of the temperature, the planetary *T*_{0}, and its amplitude *T*_{2}. Both variables evolve monotonically in the classic solution. The introduction of response time implies a coupling between *T*_{2} and the heat flux, but without any influence on *T*_{0}. Thus, the monotonical solution for *T*_{0} is kept, but with possible damped oscillations in *T*_{2} and **J**_{q}, reaching the same steady state obtained by the classical solution.

For an earth with ice caps, the classic solution can behave nonmonotonically because *T*_{0} is linked to *T*_{2} through the condition of the ice-sheet edge. Nevertheless, the nonlinearity is not enough to produce nondamped solutions. The introduction of delay time allows self-sustained oscillations beyond a certain critical value of this parameter (*τ*_{qc} ≈ 366 yr, similar to that assumed for deep oceanic circulation). However, the origin of nondamped oscillations is due to the application of both effects; the acceptance of a delay time (linking *T*_{2} and *T*_{q}) and the ice-albedo feedback (linking *T*_{0} and *T*_{2}). The existence of this periodic behavior also depends on the initial conditions assumed for the system. When the initial values of *T*_{0}, *T*_{2}, and *T*_{q} are closer to the steady-state values (*τ*_{q} = 0), the response time has to be higher in order to obtain the nondamped solution. Moreover, an increase in delay time produces a longer period for each oscillation, whereas the highest value for the planetary temperature becomes limited (a brief summary of the results found is shown in Table 2).

Bhattacharya et al. (1982) have introduced a response time for the albedo in an EBM Ghil model. These authors assumed an albedo related to past temperatures weighted through a Gaussian distribution centered at time *τ*, the delay time for the albedo. They have also incorporated *σ*, a parameter that indicates the half-width of the Gaussian distribution. With *σ* = 0, as in our case for the heat flux, they are able to obtain periodic solutions for high response times. The period found for an albedo delay time ∼1000 yr is similar to that shown in the present paper (Fig. 2).

This simple deterministic climate model can give oscillating solutions with only the internal climatic variability present by taking delay times into account in the framework of the EIT theory. In discrete one-dimensional as well as in more elaborated two- and three-dimensional models, further analysis should be made in order to observe the effects of the delay time of the dissipative fluxes in the behavior of the system. The heat of fusion, for instance, could be incorporated into the system by observing its importance in the increase of the critical response time. In this case, rough calculations indicate a *τ*_{qc} ∼450 yr, still in the range of expected values if a high decrease in oceanic circulation is accepted.

Finally, a time-dependent response time has been introduced in order to simulate strong variations in oceanic circulation, shown as a crucial element in the advent of abrupt climatic events (e.g., Stocker and Wright 1996). High variations in the temperature and the extent of the ice-sheet edge (the only possible parameters to be calculated) are in accordance with some recent studies made through complex two- and three-dimensional models (Mikolajewicz and Maier-Reimer 1990; Stocker and Wright 1996).

## Acknowledgments

This work has been partially supported by the Ministerio de Educación y Cultura Español under Contract CLI95-1867.

## REFERENCES

## Footnotes

*Corresponding author address:* Toni Pujol, Departament de Ciencies Ambientals, Universitat de Girona, Campus Montilivi, Girona, Catalonia 17071, Spain.

Email: caaps@fc.udg.es

^{1}

Furthermore, a “computational” delay time *τ*_{qcomp}, related to the grid spacing of the model should be taken into account. The total response time, therefore, would be the sum of *τ*_{q} and *τ*_{qcomp}, which can be evaluated from *τ*_{qcomp} = *l*′/*U,* where *l*′ is the cell dimension if all the cells have the same meridional length. Nevertheless, in our calculations, we only assume a delay time defined by (18) because in North’s model the variables are continuous along the latitude (*l*′ = 0) following the form of Legendre polynomials.