## 1. Introduction

Subduction is an important process for general oceanic circulation and climate variability. The dynamic nature of the subduction process, however, remains not fully understood. Often, one tends to think of the subduction of anomalies of active tracers (e.g., temperature or potential vorticity) to be the same as passive tracers. Recent theoretical and modeling studies pointed out that the subduction of active and passive tracers may have significant differences (Liu 1999; Liu and Shin 1999). This work is a further attempt to understand the observed subduction temperature anomalies in the North Pacific with the focus on the evolution of its amplitude and propagation speed.

From recent observational analyses of the North Pacific subduction temperature variability (Deser et al. 1996; Schneider et al. 1999a; Zhang and Liu 1999), one may detect that the subduction temperature anomaly propagates southward with a rapidly decaying amplitude and at a speed slower than the mean flow. Figure 1a, taken from Schneider et al. (1999a), shows the latitude–time plot of the thickness anomaly (between the 12° and 18°C isotherms) along the subduction pathway in the central North Pacific. The pathway is marked in Fig. 1b as the maximum standard deviation tongue that extends from the midlatitude southward. Even though this observed temperature anomaly is the final result of all the forcing, local and remote, wind stress and surface buoyancy fluxes, one can see that the amplitude of the variability decays rapidly southward towards 15°N (Fig. 1b), for both the warm and cold anomalies that are subducted in the 1970s and 1980s, respectively (Fig. 1a). [Schneider et al. (1999a,b) suggest that this southward weakening is part of the reason that the North Pacific subduction temperature anomaly may not significantly affect the equatorial thermocline.] The mechanism that is responsible for the amplitude reduction has not been fully studied. It is tempting to conclude that the amplitude reduction is caused simply by turbulent mixing. Here, however, as an extension of Liu (1993), we will show that the amplitude reduction is contributed to significantly by planetary wave dynamics.

Less certain is the speed of subduction. The observed temperature anomalies (Fig. 1a) appear to travel somewhat slower than the mean flow, especially north of 25°N. For example, the cold temperature anomaly patch started from 34°N in 1980 and moved to 25°N in 1990 [also seen in Figs. 10 and 16 of Deser et al. (1996)] with a southward propagation speed of about 0.3 cm s^{−1}, about half that of the mean current (about 0.7 cm s^{−1}). This slower speed of the subduction temperature anomaly has been noticed in theoretical (Liu 1999) and oceanic general circulation (Liu and Shin 1999) models and has been speculated to be related to the propagation of higher mode planetary waves. However, no clear physical explanation has been given to account for the slower speed. Why does the wave propagate always slower, rather than faster, than the ventilation flow? Here, we will show that the slower subduction speed is caused by the mean potential vorticity gradient that points eastward along the subduction pathway.

In this paper, a simple model is first used to highlight the physical mechanisms that determine the amplitude and propagation speed of a subduction wave. The theories are substantiated by observations and numerical simulations of North Pacific decadal thermocline variability. The paper proceeds as follows. The simple model is introduced in section 2. The evolution of wave amplitude and propagation speed are discussed in Sections 3 and 4, respectively. Section 5 further studies the North Pacific thermocline variability with an oceanic general circulation model. Last, conclusions and discussions are given in section 6.

## 2. Simple model

The simplest model that captures the subduction of thermocline variability is a 2-layer planetary geostrophic (PG) model (Liu 1993). The two model layers represent the stratification within the thermocline, while a rigid and flat model bottom represents the base of thermocline. This flat bottom filters out signals from the first baroclinic mode that propagates westward as a non-Doppler shift mode (Liu 1999). We therefore preferentially isolate signals of the second-baroclinic type, also known as the subduction or advective mode, that propagates downstream with the ventilation flow.

The model domain consists of the subtropical Pacific basin from 12° to 40°N. The lower layer outcrops along a latitude circle in the northern portion of the basin, while the base of the model is at *H* = 600 m. The longitudinal extent of the model domain is *L* = 8400 km, with the eastern and western boundaries specified by *x*_{e} (= 0) and *x*_{w} = −*L,* respectively.

*βHυ*

_{B}

*fw*

_{e}

*w*

_{e}is the Ekman pumping rate and

**v**

_{B}= [

**v**

_{1}

*h*+

**v**

_{2}(

*H*−

*h*)]/

*H*is the barotropic velocity that is determined by

*u*

_{B}

*υ*

_{B}

*x*

_{f}

*f*

^{2}

*w*

_{e}

*f*

*f*

^{2}

*w*

_{e}

*Hβf*

*q*

_{2}=

*f*/(

*H*−

*h*), is conserved so that

*u*

_{2},

*υ*

_{2}) are the velocities of the second layer. A single equation can be derived for the evolution of the layer interface as

*C*(

*h*) is given by

*C*

*h*

*βγh*

*H*

*h*

*f*

^{2}

*H.*

*ρ*

_{o}= 1000 kg m

^{−3}as the mean density,

*γ*=

*g*Δ

*ρ*/

*ρ*

_{o}≈ 0.02 m s

^{−2}as the reduced gravity,

*f*=

*f*

_{m}+

*βy*as the Coriolis parameter in which

*f*

_{m}= 2Ω sin(35°), and

*β*= (2Ω/

*R*) cos(35°) with

*R*as the radius of the earth.

*f*

_{o}(

*x,*

*t*), and the Ekman pumping,

*w*

_{e}(

*x,*

*f*,

*t*), respectively. This paper will focus on the thermocline anomaly produced by the surface buoyancy forcing. Therefore, the Ekman pumping remains steady according to

*w*

_{e}

*f*

*w*

_{o}

*f*

_{n}

*f*

*f*

*f*

_{s}

*f*

^{2}

_{n}

*w*

_{o}= 8.74 × 10

^{−6}m s

^{−1}(such that the maximum

*w*

_{e}is 1 × 10

^{−6}m s

^{−1}).

*f*

_{o}, the solution in the ventilated zone is simply

*h*

*H*

*f*

*f*

_{o}

*q*

_{2}=

*f*

_{o}/

*H*) within the entire ventilated zone. This is because the subduction PV is uniform along the zonal outcrop line

*f*

_{o}and in this model the depth is constant

*H*along the outcrop line. With the uniform lower-layer PV, the characteristic speeds (or group velocity) can be shown to be identical to the ventilation flow

*u*

_{B}

*C*

*h*

*υ*

_{B}

*u*

_{2}

*υ*

_{2}

**C**

_{g}

*u*

_{B}

*C*

*h*

*υ*

_{B}

## 3. Amplitude of subduction anomaly

*f*

_{o}

*x,*

*t*

*f*

_{c}

*ae*

^{−(x−xo)2/d2}

*ωt*

*a*and

*ω*are the amplitude and frequency of the perturbation outcropping latitude,

*x*

_{o}is the longitudinal location of the center of the anomaly,

*d*is the

*e*-folding distance for the zonal extent of the anomaly, and

*f*

_{c}is the mean position of the outcrop line. We now consider an example of decadal forcing of

*ω*= 2

*π*/20 yr s

^{−1}, which is similar to the North Pacific decadal variability (Fig. 1). The depth anomaly travels southward from the outcrop line,

*f*

_{c}= 36.5°N, centered at −

*x*

_{o}/

*L*= −0.3. The northern and southern boundaries are located at

*f*

_{n}= 40°N and

*f*

_{s}= 12°N, respectively. The amplitude of the anomaly is created by the meridional displacement of the outcrop line of amplitude

*a*= 1.125 deg, which corresponds to a winter sea surface temperature anomaly of approximately 1°C in the midlatitudes. The

*e*-folding distance in the zonal direction is

*d*= 84 km. Figure 2 shows several snapshots of the depth anomaly. It is seen that the depth anomaly travels along the subduction pathway, first southward and later southwestward toward the western boundary. The amplitude of the anomaly, however, seems to be reduced. The amplitude reduction can be seen more clearly in Fig. 3a, which shows the amplitude of the decadal harmonic of the depth anomaly. The most striking feature is a 40% amplitude decay downstream of the outcrop line, from 16 m to approximately 10 m at 23°N.

*T*

_{r}

*e*

^{−(x−xo)2/d2}

*ωt*

*f*

_{c}. Figure 3b shows the spatial distribution of the amplitude of the decadal harmonic of the passive tracer. In contrast to the height anomaly, the tracer maximum remains essentially unchanged downstream, as expected from (2.12).

The amplitude evolution of the depth and the passive tracer can also be seen in Figs. 4a and 4b, which show the latitude–time variation of the zonally averaged anomalies within the two bounding characteristics of subduction. The data are shown for the last 40 years of a 200-yr simulation. The basin crossing time of the anomalies is approximately 12 years. The height anomaly shows decay, while the passive tracer conserves its amplitude, consistent with Fig. 3.

The amplitude evolution in Fig. 3 is insensitive to forcing frequencies. In the extreme case of steady-state solutions, the difference of the two solutions of different outcrop lines and tracer sources gives an amplitude evolution similar to those in Fig. 3 (not shown).

*μ*= ∇ ·

**C**

_{g}is the divergence of the group velocity. In Eq. (3.4),

*μ*≠ 0, which is in sharp contrast to the passive tracer equation (2.12). Indeed, considering Eq. (2.7) for the special case of a mean zonal outcrop line (as is the case in the above example), Eq. (3.4) is identical to (2.12) except for

*μ*≠ 0 in Eq. (3.4). While the amplitude of the passive tracer is conserved downstream, the amplitude of the interface anomaly decreases (increases) for a positive (negative)

*μ.*

The group velocity, shown in (2.9), has two parts: the barotropic flow and the Rossby wave speed. Therefore, the effective damping also has two parts; *μ* = *μ*_{B} + *μ*_{R}, where *μ*_{B} = ∇ · _{B} = −*w*_{e}/*H* is the divergence of the barotropic velocity field, and *μ*_{R} = ∂_{x}*C*(*h**dC*/*d**h*_{x}*h**μ*_{B} > 0. In addition, from (2.4) we have *d*_{h}*C* < 0 in the upper thermocline (*h**H*/2). Since usually the mean depth of thermocline decreases towards the east (∂_{x}*h**μ*_{R} ≥ 0. Therefore, in a subtropical gyre, the effective damping is positive, *μ* > 0, and the thermocline anomaly decays as it propagates southward.

*f*

_{c}and therefore the mean thermocline depth in the ventilated zone is

*h*

*H*(1 −

*f*/

*f*

_{c}). This leads to ∂

_{x}

*h*

*μ*

_{R}= 0. The damping is now produced purely by the divergent barotropic flow as

*μ*=

*μ*

_{B}= −

*w*

_{e0}/

*H.*Using the Sverdrup relation (2.1), the damping timescale can be estimated as

*t*

_{D}

*H*

*w*

_{e}

*f*

*βυ*

_{B}

*R*

*υ*

_{B}

*t*

_{A}

*R*

*L*

_{y}

*L*

_{y}is the meridional width of the gyre,

*t*

_{A}=

*L*

_{y}/

*υ*

_{B}is the advection timescale in a thermocline gyre, and we have used

*β*≈

*f*/

*R.*After the gyre is completely ventilated, which has a timescale of the order of the advective timescale

*t*

_{A}, the amplitude of the thermocline anomaly is reduced to

*e*

^{−tA/tD}

*e*

^{−Ly/R}

*R*/

*L*

_{y}ranges between 2 and 3, the amplitude decays by about 30% to 40%, consistent with Fig. 3.

*a.*In light of Eq. (2.6), an anomalous change of the interface depth is given by

*f*

_{c}to the southern edge of the gyre

*f*=

*f*

_{c}−

*βL*

_{y},

*δh*decreases from

*aH*/

*f*

_{c}to

*aHf*/

*f*

^{2}

_{c}

It is important to point out that the finite changes of mean thermocline depth and Coriolis parameter are critical for the amplitude change of the subduction planetary wave. In a 2-layer quasigeostrophic (QG) thermocline model, such as that in Dewar (1989), the amplitude of a subduction planetary wave remains unchanged. Furthermore, the change of wave amplitude is independent of the variation of the PV anomaly. These two points can be demonstrated in the case of uniform mean PV.

*q*

_{2}= 0, the PV anomaly is conserved along the mean ventilation flow,

*q*

_{2}=

*f*/(

*H*−

*h*), and therefore

*q*

^{′}

_{2}

*dq*

_{2}/

*d*

*h*

*h*′ ≈

*q*

^{2}

_{2}

*h*′/

*f*. As the wave propagates,

*q*

^{′}

_{2}

*q*

_{2}=

*f*/(

*H*−

*h*

*h*′ −

*f*decreases southward in a subtropical gyre. [Note that the decay of amplitude with

*f*is consistent with the estimation in (3.6), because (3.6) was derived with uniform PV.]

For a QG model, however, *q*_{2} = −*Fϕ*_{c}, where *F* = *f*^{2}_{m}*H*_{2}*γ* and *ϕ*_{c} is the difference of the streamfunctions between layer 2 and layer 1 and is proportional to the interface deviation *h.* Therefore, *q*^{′}_{2} ≈ *ϕ*^{′}_{c}*h*′ remains unchanged during the propagation as long as the background PV is uniform, ∇*q*_{2} = 0.

## 4. Speed of thermocline anomaly

*h*′ and a passive tracer propagate at the same speed. This is a general conclusion regardless of the model and the form of PV, as long as the mean PV is uniform. Denoting the perturbation and mean fields with prime and in capitals, the linearized PV equation can always be derived (from the conservation of total PV) in the form similar to (3.7) as

*Q*would allow the perturbation PV

*q*′ to be advected by the mean flow (

*U,*

*V*), which travels at exactly the same speed as a passive tracer. It is therefore clear that the speed of the thermocline anomaly, relative to the ventilation flow, depends critically on the mean PV gradient.

As shown in Eq. (4.1), the wave propagation speed contains two parts: the mean advection [involving the first term in (4.1)] *U,* and a general-beta-induced propagation [involving the second term in (4.1)] denoted as *C*_{Q}. The latter is the speed difference between an active and a passive tracer. To estimate *C*_{Q}, however, a specific model is needed. Mathematically, a relation between the perturbation velocity (*u*′, *υ*′) and perturbation PV *q*′ needs to be prescribed before *C*_{Q} can be derived. Physically, this implies that *C*_{Q} depends not only on the mean flow and mean PV field, but also on the structure of the perturbation itself. One example is the case of the 2-layer PG model for which *C*_{Q} is estimated in the appendix.

In spite of the difficulty in estimating *C*_{Q} for the general cases, we can gain some insight by considering the basic physical argument of Rossby wave propagation. In general, with the mean PV gradient ∇*Q* pointing “north,” the general beta effect associated with ∇*Q* should propel thermocline anomalies towards the “west” (relative to the mean flow advection), like a general Rossby wave packet. In a subtropical gyre, such as the North Pacific, the mean PV gradient in the ventilated zone points predominantly southeastward [see Fig. 1b; also see Talley (1985)]. As a result, the general beta effect propels thermocline anomalies toward the northeast at a velocity *C*_{Q}, opposite to the southwestward subduction flow *U* (see Fig. 5 for a schematic figure). Thermocline anomalies therefore travel at *U* + *C*_{Q}, which is slower than the mean ventilation flow or the passive tracer. In the shadow zone, however, the mean PV gradient points northwestward, inducing a southwestward propagation *C*_{Q}, in the absence of a mean ventilation flow (*U* = 0). A thermocline anomaly that is forced, for example, by eastern boundary ventilation, will propagate southwestward with *C*_{Q}. In the pool, there is no mean PV gradient and in turn the beta-induced propagation, *C*_{Q} = 0. Therefore, both the thermocline anomaly and passive tracer travel at the same speed *U.* These three regimes of wave speed have been discussed by Liu (1999) in the context of a 2.5-layer model.

*C*

_{Q}estimate in the 2-layer model for a crude estimate of the wave propagation under realistic ocean conditions. As shown in the appendix, the wave speed along direction

*s*can be estimated as

_{n}

*q*

_{2}

*s.*To apply Eq. (4.2) to the observed subtropical North Pacific subduction temperature anomaly between the 12° and 18°C isotherms, we will choose in the 2-layer model

*H*= 600 m,

*h*

*f*= 10

^{−4}s

^{−1}. Along the subduction pathway, the mean zonal PV gradient ∂

_{x}

*q*

_{2}

_{x}

*Q*= ∂

_{x}

*q*

_{2}

*ρ*/

*ρ,*where ∂

_{x}

*Q*=

*f*Δ

*ρ*/

*ρ*Δ

*Z*∼ 2 × 10

^{−10}m

^{−1}s

^{−1}/3 × 10

^{6}m as estimated in Fig. 1b. With (4.2), these parameters lead to a northward

*C*

_{Q}∼ 0.3 cm s

^{−1}, which is about half of the mean ventilation flow speed. It is important to point out that the precise value is not important here because of the uncertainties in various approximations and parameters used in the estimation. Nevertheless, this crude estimate seems to be useful in suggesting that the general-beta-induced speed

*C*

_{Q}can be of the same order as the mean subduction flow and therefore could be responsible for a substantial speed reduction, as for the temperature anomaly in the North Pacific near 30°N (Fig. 1a).

In general, the subduction temperature anomaly should propagate slower than the mean ventilation flow in a subtropical gyre. This follows because, in the ventilated zone, the mean PV increases eastward (or, more generally, to the left of the subduction flow in the Northern Hemisphere), which leads to a poleward general-beta-induced propagation, opposite to the subduction flow. The eastward mean PV gradient can be observed in the Atlantic and Pacific (Keffer 1985; Talley 1985). This PV gradient can be understood in the ventilated thermocline model, with, say, 2.5 layers or 3.5 layers (Luyten et al. 1983). With a zonal outcrop line, the layer thickness increases to the west due to the accumulation of the downward Ekman pumping forcing. The PV on the outcrop line therefore decreases westward. This initial PV is then subducted to form the ventilated zone, resulting in a westward decrease of subduction PV in the entire ventilated zone.

We will present an example in the 2-layer model to further illustrate the difference in the speed of the subduction wave and the ventilation flow. In our 2-layer model, a nonzonal outcrop line is needed to create a PV gradient in the ventilated zone. Otherwise, as seen in (2.6), the PV along the outcrop line and, in turn, in the entire ventilated zone is constant; this leads to a total wave speed *C*_{g} the same as the subduction flow as shown in (2.7), or implied by (3.8). We will use an outcrop line that slopes northeastward in the eastern half of the basin, and northwestward in the western half of the basin. We choose this form of the outcrop to show how both a positive and negative potential vorticity gradient affects the wave speed compared to the ventilation flow.

Figure 6a shows the mean PV in the 2-layer model for our choice of the outcrop line. The potential vorticity gradient points westward west of the center PV line (of a value of about 2.55) and points eastward east of the center PV line. Figure 6b shows the difference in meridional speed between the subduction planetary wave (*C*_{gy} = *υ*_{B}) and the mean flow (*υ*_{2}), *υ*_{B} − *υ*_{2}. To the west of the centerline the wave speed *C*_{gy} (or *υ*_{B}) is more southward than the ventilation flow *υ*_{2}; to the east of the centerline the wave speed *C*_{gy} is less southward than the ventilation flow. Therefore, a subduction wave propagates southward slower (faster) than the mean flow in the eastern (western) half of the basin, where an eastward (westward) PV gradient exists.

## 5. Numerical experiments

The observed thermocline variability (Fig. 1) is subject to wind stress and surface buoyancy forcing from both local and remote regions. To better understand the evolution of temperature anomaly and its relation to the mean ventilation flow, we will perform sensitivity experiments in an oceanic general circulation model (MOM2: Pacanowski 1996). The model has a resolution of 2° × 2° in the horizontal with a total of 31 vertical layers. The model domain is the Pacific basin north of 40°S. The model seasonal climatology is spun up with the Comprehensive Ocean–Atmosphere Data Set (COADS: da Silva et al. 1994) winds, sea surface temperature (SST), and sea surface salinity (SSS). The SST and SSS forcings are imposed as a surface restoring with a restoring time of 5 days. A control simulation is performed by adding the anomalous COADS wind stress and SST from year 1945 to 1994. Choosing a band of subduction pathway similar to that in Fig. 1b, the latitudinal evolution of the subduction temperature anomaly shown in Fig. 7a (derived from White 1995) is similar to that of the thickness anomaly in Fig. 1b. The evolution of the subduction temperature anomaly in the model (Fig. 7b) bears some resemblances to the observation. Both the observations and the model show the propagation of subducted warm subsurface anomalies from the middle to low latitudes through the mid-1970s and subsequent cold anomalies through the mid-1980s. The mean PV field between 25 and 26 *σ*_{t} in the model is also shown in Fig. 8, which agrees largely with the observed PV field in Fig. 1b. Most important to our study here is the southeastward PV gradients along the subduction pathway (running roughly from 35°N, 160°E to 20°N, 180°E) in both the model (Fig. 8) and observation (Fig. 1b). There are also differences between the model simulation and observations. The magnitude of the model subduction anomaly is weaker in the subtropics, but stronger in the Tropics. The warm temperature anomaly seems to propagate faster and the cold temperature anomalies seem to penetrate less southward than in observations. The southward PV gradient in the model seems to be smaller than in observation. The PV field near the eastern boundary differs significantly from the observation, which, nevertheless, is unlikely to have much impact on our discussion of subduction. The reasons for these model–observation discrepancies are not completely clear to us. Deficiencies in the model forcing field and data could all contribute to these discrepancies. In spite of these possible model deficiencies, the reasonable resemblance between the model and observations suggests that the model is useful in simulating some major features of North Pacific decadal thermocline variability.

To focus on the subduction of a temperature anomaly that is forced by anomalous surface buoyancy forcing, we performed a sensitivity experiment (Buoyancy Run) in which the only anomalous forcing is the interannual COADS SST. For comparison with the mean ventilation flow, we also applied a passive tracer that is forced by a surface source that is identical to the anomalous restoring SST. Figure 9a shows the latitude–time plot of the temperature anomalies along the subduction pathway for the Buoyancy Run. The subduction of the temperature anomalies resembles that of the Control Run (Fig. 7b) until about 20°N. The magnitude of the subduction temperature anomalies is, however, almost 40% weaker in the sensitivity experiment. For example, the positive (negative) maximum of subduction temperature anomalies in the 1970s (1980s) are reduced from over 0.8°C (−0.8°C) in the control run to less than 0.6°C (−0.6°C) in the Buoyancy Run. South of 20°N, the Buoyancy Run produces little variability—a topic to be returned to later.

Comparing the evolution of subduction temperature anomalies in the Buoyancy Run (Fig. 9a) with the passive tracer (Fig. 9b), it is clear that the subduction temperature anomalies decay faster. For the positive anomaly that subducts in 1975, both the temperature anomaly and the passive tracer have a maximum of over 0.5 immediately after the subduction near 35°N. The temperature anomaly, however, decays rapidly to less than 0.1 before it reaches 20°N. In contrast, the 0.1 contour of the tracer extends almost to 15°N. The temperature anomaly weakens faster than the passive tracer by over 40%. Similar conclusions can be drawn for the cold anomalies that subducted in the 1980s.

The slower speed of the temperature anomaly seems to be also discernible, but less clear than the amplitude decay. The positive temperature anomaly propagates from 33°N in 1975 to 27°N in 1978, giving an estimated speed of about 0.7 cm s^{−1}.^{1} The corresponding positive passive tracer propagates from 33°N in 1975 to 25°N before 1978, giving a faster speed of 0.9 cm s^{−1}. A similar comparison can also be made for the two negative temperature anomaly pulses that start from 1983 and 1987 (Fig. 9a). To compare the propagation speeds of the temperature anomalies and the passive tracer more clearly, three characteristic lines are drawn for the temperature anomaly pulses starting from 1975, 1983, and 1987 (solid lines in Fig. 9a) by simply tracking the maximum amplitude tongue associated with each of the anomalies. The corresponding lines are also plotted for the passive tracer (dashed lines in Fig. 9b) and are superimposed on the temperature anomaly plot (dashed lines in Fig. 9a). The slower speed of the temperature anomaly can now be seen immediately. Applying (4.2) to the model simulation (in Fig. 9) would lead to a *C*_{Q} comparable to that for the observation [see discussion after (4.2)]. Considering the uncertainties of (4.2), the estimate is perhaps only useful in implying a speed reduction that is not negligible in the model (as it does for the observation), which is consistent with the discussion of Fig. 9.

In short, the major conclusions of the theoretical studies in sections 3 and 4 are supported by our model simulations in the North Pacific. Compared with the theoretical results in the last two sections, the MOM2 simulations have a much faster decay of amplitude in both the passive tracer and the temperature anomaly. To the first order, this faster decay in MOM2 is caused by model diffusion. In the real ocean, diffusive mixing is always present. Observations have shown a significant spread of passive tracer along and across isopycnal surfaces (Ledwell et al. 1993), presumably by small-scale eddies. This mixing will undoubtedly reduce the amplitude of the passive tracer and qualitatively add to the amplitude decay of the temperature anomaly, as in the MOM2 simulation. Therefore, both the mixing effect and the dynamic effect are important for the observed decay of the subduction temperature anomalies. The dynamic contribution to the total amplitude decay, as estimated before, could reach over 40%. The strength of the mixing effect in the real ocean, however, still remains to be further studied. It is likely that our MOM2 simulations overestimated the diffusive effect, particularly at the present coarse resolution.

## 6. Summary and discussions

The decadal evolution of subduction temperature anomalies in the subtropical North Pacific is studied using both a simple and a complex ocean model. It is found that the temperature anomalies evolve differently from a passive tracer in two ways. First, the amplitude of the temperature anomaly decays faster than the passive tracer by about 30%–50%. The faster decay is caused by the divergence of the group velocity of the subduction planetary wave, which is contributed to significantly by the divergent Sverdrup flow in the subtropical gyre. Second, the temperature anomaly appears to propagate slower than the passive tracer (or mean ventilation flow). The slower propagation occurs because the mean PV gradient in the ventilated zone is dominantly eastward; the associated general beta effect tends to produce a northward propagation of the temperature anomaly, as a Rossby wave response, which cancels part of the southward advection of the mean ventilation flow.

We believe that the faster decay of the temperature anomaly in the subtropics is a robust result and it has been confirmed by both theory and many sensitivity model experiments. This is also demonstrated in an accompanying study (Shin and Liu 2000) with different model configurations in both MOM2 and the Miami Isopycnal Model. At present, however, it remains to us a challenge to derive an accurate estimate of the subduction speed under realistic oceanic conditions. The idea seems straightforward: a mean eastward PV gradient would slow down the temperature anomaly relative to the mean flow. However, our attempt to estimate this speed slowdown in both observation and ocean general circulation models produces less clear results, as discussed regarding Figs. 1a, 7, and 9. Possible reasons include the oversimplification of the 2-layer model that is used to estimate the subduction speed (see appendix) and the complex surface forcing field, which may produce perturbations interacting with each other. Further studies are clearly needed for a better estimate of the speed.

Our discussions have focused on the temperature anomaly generated by surface heat flux anomaly, mostly because this anomaly can be directly compared with a passive tracer and therefore is easier to understand. It is important to recognize that thermocline temperature variability can also be generated by wind stress forcing and freshwater flux forcing. We previously speculated, through a comparison of Figs. 9a and 7b, that the wind stress variability contributes substantially to the temperature variability in the subtropics, and is responsible for most of the temperature variability in the Tropics. This speculation is further confirmed by another sensitivity experiment (Wind Run), in which only the interannual variability of COADS wind stress is applied (over the entire model domain). Comparing the latitude–time plot of the Wind Run (Fig. 10) with the corresponding plot of the control run (Fig. 7b) and the Buoyancy Run (Fig. 9a), it is clear that tropical (south of 20°N) thermocline variability is generated predominantly by the wind stress forcing. We have carried out additional sensitivity experiments in which wind anomalies are prescribed in different regions (not shown). The results show that most of the tropical variability is forced by local wind over the Tropics through baroclinic response (not through remote subduction). This point has been made clear by Schneider et al. (1999a,b). Furthermore, about half of the subduction temperature variability in the subtropics (north of about 20°N) is also generated by the wind. This may seem somewhat surprising because previous studies suggest that the subduction temperature anomaly is produced mainly by the surface buoyancy forcing (Liu 1999; Huang and Pedlosky 1999; Schneider et al. 1999a). The Ekman pumping forcing variability seems to be the most efficient in generating the first baroclinic mode, rather than the subduction planetary wave. However, it should be pointed out that these previous works have focused on the Ekman pumping effect on the permanent thermocline. Little attention has been paid to the mixed layer and its interaction with the thermocline. The wind stress forcing, in addition to its Ekman pumping effect, also distorts the mixed layer temperature and therefore can generate anomalous buoyancy flux forcing to the permanent thermocline at the bottom of the mixed layer. Another mechanism linked to the generation of the subduction planetary wave by the wind forcing is suggested by recent studies (Inui and Liu 2001, manuscript submitted to *J. Phys. Oceanogr.*, Kubokawa and Xie 2001, manuscript submitted to *J. Phys. Oceanogr.*) The subduction temperature anomaly may be caused by the interaction of the wind-generated first mode wave and the mode water.

The lack of freshwater flux forcing may also affect our model simulation. This is because part of the temperature anomaly could be salinity compensated (Miller et al. 1998). Inclusion of the salinity-compensation effect will make the temperature anomaly propagate more like a passive tracer. This effect also remains to be quantified.

We have confined our study to the subtropics and our conclusions do not directly apply to temperature anomalies that propagate in the subpolar region and toward the equatorial region. Nevertheless, the simple model results in sections 3 and 4 can be carried out similarly for the subpolar gyre. We conclude that temperature anomalies, in the subpolar gyre, should be amplified during their propagation relative to a passive tracer because of the dominant convergent flow, and in turn the group velocity. For propagation from the subtropics toward the equatorial region, the temperature anomaly may propagate more efficiently than a passive tracer because of the effect of coastal Kelvin waves (Shin and Liu 2000), opposite to that in the subtropical gyre. All of these issues suggest that the role of thermocline subduction in decadal climate variability (Zhang et al. 1998; Schneider et al. 1999b) remains to be studied carefully.

## Acknowledgments

We thank two anonymous reviewers for comments that have improved the manuscript substantially. This study is supported by a postdoctoral fellowship from University of Wisconsin—Madison (MS), and by NSF and NOAA (ZL, HY).

## REFERENCES

da Silva, A. M., C. C. Young, and S. Levitus, 1994:

*Atlas of Surface Marine Data 1994,*Vol. 2:.*Anomalies of Directly Observed Quantities,*NOAA Atlas NESDIS 7, U.S. Department Of Commerce, 416 pp.Deser, C., M. A. Alexander, and M. S. Timlin, 1996: Upper ocean thermal variations in the North Pacific during 1970–1991.

,*J. Climate***9****,**1840–1855.Dewar, W., 1989: A theory of time-dependent thermocline.

,*J. Mar. Res***47****,**1–31.Huang, R. X., and J. Pedlosky, 1999: Climate variability inferred from a layered model of the ventilated thermocline.

,*J. Phys. Oceanogr***29****,**779–790.Keffer, T., 1985: The ventilation of the world's oceans: Maps of the potential vorticity field.

,*J. Phys. Oceanogr***15****,**509–523.Ledwell, J., A. J. Watson, and C. S. Law, 1993: Evidence for slow mixing across the pycnocline from an open-ocean tracer release experiment.

,*Nature***364****,**701–703.Levitus, S., 1982:

*Climatological Atlas of the World Ocean*. NOAA Prof. Paper No. 13, U.S. Govt. Printing Office, Washington, DC, 173 pp.Liu, Z., 1993: Thermocline forced by varying wind. Part II: Annual and decadal Ekman pumping.

,*J. Phys. Oceanogr***23****,**2523–2541.——,. 1999: Forced planetary wave response in a thermocline gyre.

,*J. Phys. Oceanogr***29****,**1036–1055.——, and Pedlosky, J., 1994: Thermocline forced by annual and decadal surface temperature variation.

,*J. Phys. Oceanogr***24****,**587–608.——, and Shin, S., 1999: On thermocline ventilation of active and passive tracers.

,*Geophys. Res. Lett***26****,**357–360.Luyten, J., J. Pedlosky, and H. Stommel, 1983: The ventilated thermocline.

,*J. Phys. Oceanogr***13****,**292–309.Miller, A., D. Cayan, and W. White, 1998: A westward-intensified decadal change in the North Pacific thermocline and gyre-scale circulation.

,*J. Climate***11****,**3112–3127.Pacanowski, R. C., 1996: MOM2 documentation. GFDL Ocean Tech. Rep. 3.2, 329 pp.

Schneider, N., A. J. Miller, M. A. Alexander, and C. Deser, 1999a: Subduction of decadal North Pacific temperature anomalies: Observations and dynamics.

,*J. Phys. Oceanogr***29****,**1056–1070.——, Venzke, S., A. Miller, D. Pierce, T. Barnett, C. Deser, and M. Latif, 1999b: Pacific thermocline bridge revisited.

,*Geophy. Res. Lett***26****,**1329–1332.Shin, S-I., and Z. Liu, 2000: On the response of the equatorial thermocline to extratropical buoyancy forcing.

,*J. Phys. Oceanogr***30****,**2883–2905.Talley, L., 1985: Ventilation of the subtropical North Pacific: The shallow salinity minimum.

,*J. Phys. Oceanogr***15****,**633–649.White, W. B., 1995: Design of a global observing system for gyre-scale upper ocean temperature variability.

*Progress in Oceanography,*Vol. 36, Pergamon, 169–217.Zhang, R. H., and Z. Liu, 1999: Decadal thermocline variability in the North Pacific: Two pathways around the subtropical gyre.

,*J. Climate***12****,**3273–3296.——, Rothstein, L. M., and A. J. Busalacchi, 1998: Origin of upper ocean warming and El Niño change on decadal scales in the tropical Pacific Ocean.

,*Nature***391****,**879–883.

## APPENDIX

### Estimate the General-Beta-Induced Wave Speed *C*_{Q}

*C*

_{Q}can be estimated as follows. First, one can show a relation between the second-layer velocity and barotropic velocity:

*u*

_{T},

*υ*

_{T}) = (

*u*

_{1}−

*u*

_{2},

*υ*

_{1}−

*υ*

_{2}) = −(

*γ*/

*f*)(−∂

_{y}

*h,*∂

_{x}

*h*) is the thermal wind. For a baroclinic disturbance, the perturbation velocity is due to the thermal wind advection and can be written as

*s, n*) along the mean streamline, where

*s*points downstream. Assuming that the curvature of the mean streamline is not very large, the perturbation PV equation (3.7) can be written as

_{t}

*U*

_{2}

_{s}

*q*

^{′}

_{2}

*υ*

^{′}

_{2}

_{n}

*q*

_{2}

*q*

^{′}

_{2}

*h*′

*f*/(

*H*−

*h*

^{2}, we can rewrite the perturbation PV equation as

*υ*

^{′}

_{2}

_{n}

*q*

_{2}therefore consists of two terms, the second and last terms in (A.2), which are associated with the advection on the mean PV gradient by the perturbation thermal wind and mean thermal wind, respectively; the former contributes to the propagation speed of

*h*′, while the latter contributes to the amplitude change of

*h*′, as is obvious in the form of (A.3). The general-beta-induced propagation speed in the direction of the mean flow is contributed by the second term in (A.3) and therefore can be estimated, in the WKB sense, as

*C*

_{Q}≈ −∂

_{n}

*q*

_{2}

*γ*

*h*

*H*−

*h*

^{2}/

*Hf*

^{2}.

Snapshots of the 2-layer model interface anomaly forced by a perturbation outcrop line at 20-yr period [Eq. (3.1)]. Other forcing parameters are described in the text. (a) *t* = *t*_{o}, (b) *t* = *t*_{o} + 5 yr, (c) *t* = *t*_{o} + 10 yr (contour interval is 2 m; negative contours are dashed)

Citation: Journal of Physical Oceanography 31, 7; 10.1175/1520-0485(2001)031<1733:EOSPWW>2.0.CO;2

Snapshots of the 2-layer model interface anomaly forced by a perturbation outcrop line at 20-yr period [Eq. (3.1)]. Other forcing parameters are described in the text. (a) *t* = *t*_{o}, (b) *t* = *t*_{o} + 5 yr, (c) *t* = *t*_{o} + 10 yr (contour interval is 2 m; negative contours are dashed)

Citation: Journal of Physical Oceanography 31, 7; 10.1175/1520-0485(2001)031<1733:EOSPWW>2.0.CO;2

Snapshots of the 2-layer model interface anomaly forced by a perturbation outcrop line at 20-yr period [Eq. (3.1)]. Other forcing parameters are described in the text. (a) *t* = *t*_{o}, (b) *t* = *t*_{o} + 5 yr, (c) *t* = *t*_{o} + 10 yr (contour interval is 2 m; negative contours are dashed)

Citation: Journal of Physical Oceanography 31, 7; 10.1175/1520-0485(2001)031<1733:EOSPWW>2.0.CO;2

(a) Amplitude of the decadal interface (*h*) variability in Fig. 2. (b) Amplitude of decadal passive tracer variability forced by (3.2). The *h* amplitude is normalized by its maximum (16 m)

Citation: Journal of Physical Oceanography 31, 7; 10.1175/1520-0485(2001)031<1733:EOSPWW>2.0.CO;2

(a) Amplitude of the decadal interface (*h*) variability in Fig. 2. (b) Amplitude of decadal passive tracer variability forced by (3.2). The *h* amplitude is normalized by its maximum (16 m)

Citation: Journal of Physical Oceanography 31, 7; 10.1175/1520-0485(2001)031<1733:EOSPWW>2.0.CO;2

(a) Amplitude of the decadal interface (*h*) variability in Fig. 2. (b) Amplitude of decadal passive tracer variability forced by (3.2). The *h* amplitude is normalized by its maximum (16 m)

Citation: Journal of Physical Oceanography 31, 7; 10.1175/1520-0485(2001)031<1733:EOSPWW>2.0.CO;2

Latitude–time plot of zonally averaged decadal (a) *h* and (b) passive tracer anomaly in the 2-layer model. The *h* anomaly is normalized by its maximum of 16 m. (Negative contours are dashed.)

Citation: Journal of Physical Oceanography 31, 7; 10.1175/1520-0485(2001)031<1733:EOSPWW>2.0.CO;2

Latitude–time plot of zonally averaged decadal (a) *h* and (b) passive tracer anomaly in the 2-layer model. The *h* anomaly is normalized by its maximum of 16 m. (Negative contours are dashed.)

Citation: Journal of Physical Oceanography 31, 7; 10.1175/1520-0485(2001)031<1733:EOSPWW>2.0.CO;2

Latitude–time plot of zonally averaged decadal (a) *h* and (b) passive tracer anomaly in the 2-layer model. The *h* anomaly is normalized by its maximum of 16 m. (Negative contours are dashed.)

Citation: Journal of Physical Oceanography 31, 7; 10.1175/1520-0485(2001)031<1733:EOSPWW>2.0.CO;2

The schematic figure illustrating the propagation speed of a thermocline anomaly. The thin solid lines are mean PV contours; *U* and *C*_{Q} are the mean advection and general-beta-induced wave propagation velocity: P.Z., V.Z., and S.Z. denote the pool zone, the ventilated zone and the shallow zone, respectively

Citation: Journal of Physical Oceanography 31, 7; 10.1175/1520-0485(2001)031<1733:EOSPWW>2.0.CO;2

The schematic figure illustrating the propagation speed of a thermocline anomaly. The thin solid lines are mean PV contours; *U* and *C*_{Q} are the mean advection and general-beta-induced wave propagation velocity: P.Z., V.Z., and S.Z. denote the pool zone, the ventilated zone and the shallow zone, respectively

Citation: Journal of Physical Oceanography 31, 7; 10.1175/1520-0485(2001)031<1733:EOSPWW>2.0.CO;2

The schematic figure illustrating the propagation speed of a thermocline anomaly. The thin solid lines are mean PV contours; *U* and *C*_{Q} are the mean advection and general-beta-induced wave propagation velocity: P.Z., V.Z., and S.Z. denote the pool zone, the ventilated zone and the shallow zone, respectively

Citation: Journal of Physical Oceanography 31, 7; 10.1175/1520-0485(2001)031<1733:EOSPWW>2.0.CO;2

(a) Mean PV in a 2-layer model thermocline forced under a nonzonal outcrop line. The outcrop line is marked as heavy lines. (b) The normalized meridional velocity difference of barotropic and ventilation flow (*υ*_{B} − *υ*_{2})

Citation: Journal of Physical Oceanography 31, 7; 10.1175/1520-0485(2001)031<1733:EOSPWW>2.0.CO;2

(a) Mean PV in a 2-layer model thermocline forced under a nonzonal outcrop line. The outcrop line is marked as heavy lines. (b) The normalized meridional velocity difference of barotropic and ventilation flow (*υ*_{B} − *υ*_{2})

Citation: Journal of Physical Oceanography 31, 7; 10.1175/1520-0485(2001)031<1733:EOSPWW>2.0.CO;2

(a) Mean PV in a 2-layer model thermocline forced under a nonzonal outcrop line. The outcrop line is marked as heavy lines. (b) The normalized meridional velocity difference of barotropic and ventilation flow (*υ*_{B} − *υ*_{2})

Citation: Journal of Physical Oceanography 31, 7; 10.1175/1520-0485(2001)031<1733:EOSPWW>2.0.CO;2

Latitude–time plots of annual mean temperature anomalies from (a) the XBT data of White (1995) and (b) the control experiment. The seasonal climatology from 1970 to 1994 is subtracted to derive the temperature anomalies. The temperature anomaly is averaged zonally within the subduction band similar to that marked in Fig. 1a. The contour interval is 0.2°C and the heavy shading and light shading represent regions larger than 0.2°C and smaller than −0.2°C, respectively. (Negative contours are dashed.)

Citation: Journal of Physical Oceanography 31, 7; 10.1175/1520-0485(2001)031<1733:EOSPWW>2.0.CO;2

Latitude–time plots of annual mean temperature anomalies from (a) the XBT data of White (1995) and (b) the control experiment. The seasonal climatology from 1970 to 1994 is subtracted to derive the temperature anomalies. The temperature anomaly is averaged zonally within the subduction band similar to that marked in Fig. 1a. The contour interval is 0.2°C and the heavy shading and light shading represent regions larger than 0.2°C and smaller than −0.2°C, respectively. (Negative contours are dashed.)

Citation: Journal of Physical Oceanography 31, 7; 10.1175/1520-0485(2001)031<1733:EOSPWW>2.0.CO;2

Latitude–time plots of annual mean temperature anomalies from (a) the XBT data of White (1995) and (b) the control experiment. The seasonal climatology from 1970 to 1994 is subtracted to derive the temperature anomalies. The temperature anomaly is averaged zonally within the subduction band similar to that marked in Fig. 1a. The contour interval is 0.2°C and the heavy shading and light shading represent regions larger than 0.2°C and smaller than −0.2°C, respectively. (Negative contours are dashed.)

Citation: Journal of Physical Oceanography 31, 7; 10.1175/1520-0485(2001)031<1733:EOSPWW>2.0.CO;2

The mean model PV field between 25–26 *σ*_{t} surfaces similar to Fig. 1b. The contour interval is 10^{−10} m^{−1} s^{−1}. The PV is derived from the run that is forced by the climatological COADS surface forcing

Citation: Journal of Physical Oceanography 31, 7; 10.1175/1520-0485(2001)031<1733:EOSPWW>2.0.CO;2

The mean model PV field between 25–26 *σ*_{t} surfaces similar to Fig. 1b. The contour interval is 10^{−10} m^{−1} s^{−1}. The PV is derived from the run that is forced by the climatological COADS surface forcing

Citation: Journal of Physical Oceanography 31, 7; 10.1175/1520-0485(2001)031<1733:EOSPWW>2.0.CO;2

The mean model PV field between 25–26 *σ*_{t} surfaces similar to Fig. 1b. The contour interval is 10^{−10} m^{−1} s^{−1}. The PV is derived from the run that is forced by the climatological COADS surface forcing

Citation: Journal of Physical Oceanography 31, 7; 10.1175/1520-0485(2001)031<1733:EOSPWW>2.0.CO;2

Latitudinal evolution of subduction (a) temperature and (b) passive tracer anomalies similar to Fig. 7, but for the Buoyancy Run. The contour interval is now 0.1

Citation: Journal of Physical Oceanography 31, 7; 10.1175/1520-0485(2001)031<1733:EOSPWW>2.0.CO;2

Latitudinal evolution of subduction (a) temperature and (b) passive tracer anomalies similar to Fig. 7, but for the Buoyancy Run. The contour interval is now 0.1

Citation: Journal of Physical Oceanography 31, 7; 10.1175/1520-0485(2001)031<1733:EOSPWW>2.0.CO;2

Latitudinal evolution of subduction (a) temperature and (b) passive tracer anomalies similar to Fig. 7, but for the Buoyancy Run. The contour interval is now 0.1

Citation: Journal of Physical Oceanography 31, 7; 10.1175/1520-0485(2001)031<1733:EOSPWW>2.0.CO;2

Latitudinal evolution of subduction temperature anomalies similar to Fig. 9a, but for the Wind Run

Citation: Journal of Physical Oceanography 31, 7; 10.1175/1520-0485(2001)031<1733:EOSPWW>2.0.CO;2

Latitudinal evolution of subduction temperature anomalies similar to Fig. 9a, but for the Wind Run

Citation: Journal of Physical Oceanography 31, 7; 10.1175/1520-0485(2001)031<1733:EOSPWW>2.0.CO;2

Latitudinal evolution of subduction temperature anomalies similar to Fig. 9a, but for the Wind Run

Citation: Journal of Physical Oceanography 31, 7; 10.1175/1520-0485(2001)031<1733:EOSPWW>2.0.CO;2

^{1}

The faster speed of model subduction anomaly in the Buoyancy Run than in the observation (Fig. 1a or 7a) could be due to model deficiency. It could also be caused by locally forced responses in the observed thermocline. Indeed, when the full forcing is applied, the subduction propagation in the Control Run (Fig. 7b) does not shown a systematic bias toward a higher speed than in the observation (Fig. 7a).