Climatic variability of the pycnocline in the eastern tropical and North Pacific has oceanographic and ecological implications. Gridded monthly profiles of temperature and salinity from the Simple Ocean Data Assimilation (SODA) reanalysis, 1958–2008, were used to derive estimates of four variables related to the density structure of the upper-ocean water column: surface temperature, pycnocline depth, mixed layer depth, and stratification (potential energy anomaly). The pycnocline is primarily a thermal gradient in this region, except in subarctic waters at the northern extreme of the study area, where salinity becomes more important than temperature in determining stratification. Spatial patterns of mean and standard deviation of the four pycnocline variables are presented. Partitioning of variance between seasonal and interannual scales shows the predominance of interannual variability in the tropics and seasonal variability at higher latitudes. Low-frequency variations (trends) in the pycnocline variables were derived by state-space analysis of time series averaged in 5° squares. Regionally coherent trends were either monotonic over 50 years or had decadal-scale changes in sign (±5–10-m depth, ±5%–10% of stratification). For example, in the eastern equatorial Pacific, the pycnocline shoaled by 10 m and weakened by 5% over the 50 years, while in the California Current the pycnocline deepened by ~5 m but showed little net change in stratification, which weakened by 5% to the mid-1970s, strengthened by 8% to the mid-1990s, and then weakened by 4% to 2008. These observed changes in the pycnocline, and future changes resulting from global climate change, may have important biological and ecosystem effects.
The pycnocline is a layer of water in which density changes most rapidly with depth (Sprintall and Cronin 2001). The pycnocline is usually associated with a temperature gradient or thermocline, because density variations with depth are determined primarily by temperature, except at high latitudes where salinity is more important. The pycnocline is a physical gradient that affects buoyancy, heat budgets, currents, and transport. The pycnocline is also an ecological boundary both because it may include a physiological temperature limit and because it often corresponds to gradients in nutrients, oxygen, or other limiting factors. This “ubiquitous horizontal front” is “globally associated with the change from epipelagic to deeper ecosystems, so it is the most significant feature in the three-dimensional ecological geography of the sea” (Longhurst 2007). Longhurst goes on to review how spatial and temporal variability of the pycnocline determines global patterns of primary productivity and ecosystem structure, resulting in a set of biogeochemical provinces with distinct oceanographic and ecological characteristics.
Annual variability of thermocline depth and stratification was first invoked as a determinant of the timing of spring phytoplankton blooms by Sverdrup (1953). The critical depth hypothesis has been a cornerstone of pelagic ecology since then, although it has recently been reexamined (Behrenfeld 2010). Interannual variations of the thermocline in the equatorial Pacific are fundamental to El Niño–Southern Oscillation (ENSO) dynamics (Wyrtki 1985) and force variations at higher latitudes (Miller et al. 1997). On longer time scales, McGowan et al. (2003) linked decreasing zooplankton volume in the California Current since 1951 with increasing thermocline stratification, although Lavaniegos and Ohman (2007) suggested that this was a change in the composition (decline in salps and doliolids) rather than the total biomass of the zooplankton community.
The density structure of the upper-ocean water column is determined by a number of processes including fluxes of heat and freshwater at the air–sea interface, diffusive and turbulent mixing, vertical advection of water by convection and by upwelling and downwelling processes, and horizontal movement of water or water properties by wind-driven and general circulation and by Rossby, Kelvin, and higher-frequency internal waves. The depth of the mixed layer is more closely linked to local processes, such as wind mixing, insolation, and rainfall, than is the depth of the pycnocline, which is also controlled by baroclinicity associated with the ocean circulation (Longhurst 2007). The upper-ocean pycnocline considered here (defined below as the maximum density gradient) corresponds to the seasonal pycnocline underlying the surface mixed layer (Sprintall and Cronin 2001) rather than to the deeper permanent pycnocline that separates stratified upper-ocean waters from dense abyssal waters (e.g., Gnanadesikan 1999).
These processes controlling the structure of the pycnocline are neither constant in time nor uniform in space. In this paper, we analyze spatial and temporal patterns of variability of sea surface temperature, mixed layer depth, pycnocline depth, and stratification in a region of the eastern and central tropical and northern Pacific (Fig. 1). This region covers ecosystems and fisheries monitored and/or managed by the U.S. National Marine Fisheries Service (http://swfsc.noaa.gov/) and the California Cooperative Oceanic Fisheries Investigations (http://www.calcofi.org/), among other regulatory and academic entities. Surface temperature is included in the analysis primarily because it is a variable that is often used to depict spatial pattern and temporal change in the marine environment (e.g., Cane et al. 1997; Enfield and Mestas-Nuñez 2000).
An initial analysis will separate variance into seasonal and interannual (>1 yr) components. The focus will thereafter be on longer-term variability, which is usually categorized as decadal (or interdecadal or multidecadal) variability and as linear trends. Here, following the state-space analytical method described below, we will refer to the longer-term mode of variability in a time series as the “trend component.” Trends represent a baseline about which seasonal and interannual cyclic variations occur; trends may be time-varying or constant. For example, the gradual and continually changing background state for cyclic ENSO variations of surface temperature (Philander and Federov 2003) is a trend in a state-space analysis of an eastern equatorial Pacific SST time series (Mendelssohn et al. 2005). We will ask whether there are regionally coherent trends in surface temperature and pycnocline variables and how these trends are related among variables. The seasonal and interannual cyclic components of the state-space results will be the subject of a future paper. Organisms in both terrestrial and marine environments have evolved developmental, reproductive, and migratory strategies in adaptation to such periodic variability (Overland et al. 2010). The unpredictability of longer-term variations, as a baseline for seasonal and interannual cyclic variations, means that these changes may have effects disproportionate to their relative amplitude. Rice (2001) suggested that the effect of environmental variations on ecosystem dynamics depends on their temporal scale relative to the generation time of important species such as top predators. Our focus in this paper is to describe these longer-term variations and to discuss implications for ecosystem productivity and processes.
2. Data and methods
Monthly three-dimensional fields of temperature and salinity were acquired from the Simple Ocean Data Assimilation (SODA) 2.1.6 dataset. SODA is a reanalysis of ocean climate using a global ocean circulation model forced by observed surface wind stresses and constrained by constant assimilation of observed temperatures, salinities, and altimetry using an optimal data assimilation technique (Carton and Giese 2008). Carton et al. (2000) conducted a comparison study using independent observations to examine the accuracy of the SODA analysis. Numerous published studies have used SODA output to investigate the variability in physical ocean conditions in various ocean basins (Shi et al. 2007; Ashok et al. 2004; Moon et al. 2004). SODA 2.1.6 is the most recent “released” version. The model output is monthly averaged, covers 1958–2008, and is mapped onto a uniform 0.5° × 0.5° × 40-level grid. The reanalysis provides three output variables (temperature, salinity, velocity); temperature and salinity are the assimilated variables and thus the output values are well constrained by observations. To estimate upper-ocean pycnocline variability in the eastern tropical and North Pacific area described in the introduction, SODA model output were selected from 29.75°S to 49.75°N latitude and from 179.75° to 70.25°W longitude (0.5° resolution) at the upper 23 depth levels of 5.01, 15.07, 25.28, 35.76, 46.61, 57.98, 70.02, 82.92, 96.92, 112.32, 129.49, 148.96, 171.4, 197.79, 229.48, 268.46, 317.65, 381.39, 465.91, 579.31, 729.35, 918.37, and 1139.15 m. SODA data were obtained from the National Oceanic and Atmospheric Administration (NOAA)–National Marine Fisheries Service–Southwest Fisheries Science Center Environmental Research Division Data Access Program (ERDDAP; http://coastwatch.pfeg.noaa.gov/erddap).
At each monthly latitude/longitude grid point, the temperature and salinity at depth were used to calculate potential density. Before extracting descriptive variables for the pycnocline (mixed layer depth, pycnocline depth, stratification), potential density profiles were interpolated at 1-m intervals by cubic spline interpolation between values at 24 depths, after setting the surface values equal to the 5.01-m values. Since these profiles were relatively smooth with no noise due to sampling error, pycnocline depths were estimated simply as the depth z of the maximum difference in density between depths (z − 5 m) and (z + 5 m); see Fiedler (2010) for a discussion of problems that can arise when estimating thermocline variables from observed profile data.
Mixed layer depth was estimated as the depth at which potential density has increased relative to the surface value by an amount equivalent to a 0.5°C decrease in temperature. While Kara et al. (2000) recommended ΔT = −0.8°C as a measure of mixed (isothermal) layer depth, the ΔT = −0.5°C criterion used here is equivalent to changes in potential density, Δσθ, of +0.08 to +0.16 kg m−3, which are within the range of values used by other researchers to define mixed layer depth (de Boyer Montégut et al. 2004).
Stratification was quantified as potential energy anomaly (ϕ, J m−3), which is equal to the amount of work per unit volume required to redistribute the mass in a complete mixing of a water column to a specified depth (Simpson 1981):
where ρ is the potential density at depth z. The specified depth h for this study is 300 m, to cover the pycnocline at all locations and times of year within the study area. The relative importance of temperature and salinity in determining stratification was assessed by calculating potential energy anomaly at constant (=average) salinity ϕTP, constant temperature ϕSP, and constant temperature and salinity ϕP. Thus, stratification due to temperature ϕT = ϕTP − ϕP, and stratification due to salinity ϕS = ϕSP − ϕP. The value of ϕ due to depth variations of temperature, salinity, and pressure ranged from 631 to 1855 J m−3, while the value of ϕP, which is due only to the compressibility of seawater of average temperature and salinity, was relatively constant and ranged from 317 to 341 J m−3.
The 51-yr time series at each 0.5° grid point were analyzed first to illustrate spatial patterns of the mean and of variability at seasonal and interannual time scales. In this paper, variability is quantified as standard deviation. Both standard deviation and variance depend on the shape of the frequency distribution of the variable, including its spread or range. However, recall from elementary statistics that the standard deviation represents the dispersion or range about the mean of a variable, while variance is the standard deviation squared. Thus, total variability is the standard deviation of the monthly time series. The variance of each time series was partitioned into interannual (variance of monthly anomalies) and seasonal (total variance minus variance of monthly anomalies). Spatial patterns of these temporal variabilities (standard deviation) are displayed as contour maps.
While this crude partitioning of variance is useful for describing the spatial patterns and relative importance of interannual and seasonal variability in the large geographic area covered in this study, it does not provide a description of the temporal patterns of variability. The SODA monthly time series were therefore analyzed using a state-space decomposition that has been applied previously to examine long-term changes in the mean and seasonal components in many climate time series (Schwing and Mendelssohn 1997; Mendelssohn and Schwing 2002; Mendelssohn et al. 2003, 2005, and references therein). This method separates possibly nonstationary dynamics of the data on different time scales in an objective manner. State-space models assume that each observation y(t), in a time series t = 1, …, τ, is the sum of four components,
where, at time t, T(t) is the unobserved time-dependent mean level (nonlinear trend), S(t) is the unobserved seasonal component (zero mean, nonstationary, and nondeterministic), I(t) is the unobserved cycle component (containing any stationary autocorrelated or cyclic parts of the data), and e(t) is the stationary uncorrelated component or “observational error.” Piecewise continuous smoothing splines are used to estimate the unobserved components [for a good introduction to state-space models, see Durand and Mendelssohn (1998) and Commandeur et al. (2011); for a more in-depth discussion, see Commandeur and Koopman (2007) and Durbin and Koopman (2001)].
Advantages of the state-space approach are that 1) the trend component is nonparametric and can capture complex dynamics; 2) the cyclic component, as shown below, can be deterministic but generally is a stationary, stochastic cycle with time-varying phase and amplitude; and 3) the seasonal component has a constant period but may vary in amplitude. In more detail, the trend term can be viewed as an unknown function of time, and is parameterized as
where ∇ is the lag difference operator; k is the degree of the lag, which is equal to 1 for all our analyses, and N(0, ) denotes a random variable that is normally distributed with mean 0 (zero) and variance , which is estimated and controls the smoothness of the estimated trend. The trend gives a nonparametric estimate of the change of the level (mean) of the series with time.
We constrain the running sum of the seasonal component S as follows:
where is estimated and controls the smoothness of the estimated seasonal component, and s = 12 for monthly data.
The state-space specification of the irregular cyclic term is
where and are the states, λc is the frequency, in radians, in the range 0 < λc < π, and are two mutually uncorrelated white noise disturbances with zero means and common variance , and ρ is a damping factor. The damping factor ρ accounts for the time over which a higher amplitude event (consider this to be a “shock” to the series) in the stochastic cycle will contribute to subsequent cycles. A stochastic cycle has changing amplitude and phase, and becomes a first-order autoregression if λc is 0 or π. Moreover, it can be shown that as ρ→1, then and the stochastic cycle reduces to the stationary deterministic cycle:
The observation errors are assumed to be zero mean, independent, and identically distributed as
The state-space model is fit to each time series with a Kalman filter/smoother (cf. Commandeur et al. 2011). Maximum likelihood estimates of the unknown parameters in the model are obtained by maximizing the log likelihood for a given ψ, computed using the output of the Kalman filter. The final estimates of each component (trend, seasonal, cyclic) are obtained from the Kalman-smoothed values calculated at the maximum likelihood estimate of the parameters. Checks are made that the model assumptions are met by examining properties of the predicted and smoothed residuals. More information on the methods can be found in volume 41 of the Journal of Statistical Software, a special volume devoted to state-space models and available online (at http://www.jstatsoft.org/v41).
Time series were averaged in 136 five-degree squares, selected to cover the study area with more complete coverage in the more variable coastal regions and less complete coverage in the more homogeneous subtropical gyres. (The spatial coverage by selected squares can be seen in the lower right panel of Fig. 5, for example, where the nonwhite squares mark the locations of the analyzed time series.) This selection of 5° squares was necessary to reduce computation time—as mentioned above, an actual model (plus model checking) is estimated by maximum likelihood at each location for each variable. Random locations among the analyzed squares were selected for detailed analysis in order to determine the appropriate model for these data. A variety of models were fit, and compared using both residual analysis and model fit statistics. The model chosen for all locations included a nonparametric trend, a stochastic seasonal component (a nonstationary cycle with fixed annual frequency), and one or two stochastic cycles, whose frequencies were estimated by maximum likelihood in the modeling procedure (the mixed layer depth time series needed only one stochastic cycle, whereas the other variables were better fit with two stochastic cycles). This model did not always effectively account for variability at different time scales in each of the time series; for example, cyclic variations of a few years that were not fit by the cyclic term(s) of the model could show up in the trend term, as will be seen in a few of the trend components below. Interannual variability with multiple recurrences in the time series is captured by the cycle component(s). Longer-term variability, relative to the length of the time series, appears in the time-varying mean or trend component, which will be the focus of the state-space decomposition results.
Trend components were objectively clustered, rather than grouped by predetermined oceanographic or ecological regions (e.g., Fig. 1), in order to assess regional coherence or lack of coherence in the data. Spatial clusters for the data were calculated using the “FindClusters” command in the Mathematica software. This is a k-means type clustering algorithm that uses the median rather than the mean in computing the distances, providing clusters that are more resistant to outliers in the data. The method requires specification of the number of clusters. In trials with three, four, five, and seven clusters, five clusters were found to give both coherent and oceanographically meaningful spatial patterns for the study area. Fewer clusters lumped regions or parts of regions that are known to be fundamentally different (e.g., eastern boundary currents and subtropical gyres), while more splitting gave more spatial intermingling of squares assigned to clusters as differences in trend shapes became less distinct. Trends were normalized (mean = 0 and standard deviation = 1) prior to cluster analysis, so that each cluster contained trends of similar shape (direction of trend and timing of changes in trend), irrespective of the mean and magnitude. Illustrated trends (see Figs. 5–8) are standardized (mean = 0), but not normalized to show differences in magnitude. Time series with a constant trend (fixed mean) were excluded from clustering and were not considered in interpreting the results of this analysis.
3. Mean and variability
a. Surface temperature
Surface temperature shows expected patterns of mean and variability (Fig. 2a). The eastern Pacific warm pool, centered off southwestern Mexico, is north of the equatorial cold tongue, whereas the warm waters straddling the equator in the central Pacific correspond to the eastern edge of the western Pacific warm pool. The local temperature minimum southeast of the center of the eastern Pacific warm pool is the Costa Rica Dome. Surface temperature decreases poleward and in the eastern boundary currents. Surface temperature variability is highest in the North Pacific subarctic–subtropical transition zone (Roden 1971), the equatorial cold tongue and the Peru Current, in the subtropical South Pacific at the southern edge of the study area, and in the vicinity of southern Baja California. Interannual variability exceeds seasonal variability in the equatorial Pacific (Fig. 2a, right panel), which is known to be the center of action for ENSO dynamics. In the east, this occurs at the northern edge of the equatorial cold tongue and in the vicinity of the Costa Rica Dome. Toward the central equatorial Pacific, interannual variability becomes relatively more important as the thermocline deepens and seasonality of the cold tongue decreases.
b. Pycnocline depth
The pycnocline is generally more shallow toward the coast, associated with the eastern boundary current circulation (Fig. 2b). Zonal ridges along the equator and ~10°N are associated with equatorial surface currents (Fiedler and Talley 2006; Kessler 2006). The 10°N countercurrent ridge peaks at the Costa Rica Dome (Fiedler 2002a). The depressions in the pycnocline poleward of the equatorial current system (~20°S–20°N) correspond to the subtropical convergences of wind-driven surface currents (Tomczak 2001). Variability of pycnocline depth is greatest at these subtropical locations and least in the equatorial region and in the eastern tropical Pacific, where the permanent thermocline is strong and shallow and a seasonal thermocline does not develop (Fiedler and Talley 2006). Unlike surface temperature, interannual variability exceeds seasonal variability throughout the study area, except in the North Pacific transition zone and in the southeastern subtropical gyre. The zonal bands in the interannual fraction of variability between 20°N and 20°S are caused by seasonal variability of the equatorial circulation; the North Equatorial Countercurrent (NECC) weakens in boreal spring and wind-driven equatorial upwelling strengthens later in the year (Kessler 2006).
c. Mixed layer depth
Spatial patterns of the mean and variability of mixed layer depth are very similar to those for pycnocline depth (Fig. 2c). By definition, mixed layer depth is less than pycnocline depth; as a result, total variability is less. Seasonal variability of mixed layer depth is greater than interannual variability north of 30°N and south of 20°S. Overall, seasonality is relatively greater for mixed layer depth than it is for pycnocline depth because local processes with seasonal variations—heat exchange and wind mixing—are more important in determining mixed layer depth.
The mean pattern of stratification in the surface layer (0–300 m) is similar to the pattern of sea surface temperature: stratification is higher in warmer tropical and equatorial waters than in cooler subtropical and subarctic waters (Fig. 2d). Temperature dominates stratification of the 0–300-m surface layer except in subarctic waters at the northern extreme of the study area (Fig. 3, bottom). High-salinity Subtropical Surface Water and, closer to the equator at ~100-m depth, Subtropical Underwater (Fiedler and Talley 2006) results in a pycnocline that is weakened by a negative salinity gradient below the warm surface layer (Fig. 3, bottom, T > −S).
Like surface temperature, and unlike pycnocline and mixed layer depth, variability of stratification tends to be higher where the mean is lower (Fig. 2d). Variability is greatest in tropical surface waters along 10°N, where seasonal rainfall in the intertropical convergence zone is high (Amador et al. 2006); note the local minimum in the interannual fraction of variability along the ITCZ. Variability is also high along the equator and in the subtropical gyre northwest of Hawaii, but the interannual fraction of variability differs in these two areas. Interannual variability exceeds seasonal variability in much of the region, especially within a narrow band along the equator. Seasonal variability tends to dominate in the North Pacific transition zone and the South Pacific subtropical gyre.
The concept of time-varying or constant trends representing a baseline about which seasonal and interannual variations occur is illustrated in Fig. 4. For example, in the eastern equatorial Pacific (Fig. 4, top), surface temperature variability during 1981–2000 is primarily cyclic, with ±2°–3°C warmings and coolings associated with the 1982–85, 1986–89, and 1997–99 ENSO (El Niño and La Niña) events, along with a ±0.8°C seasonal cycle and a small trend component comprising a warming of 0.1°C. In the California Current (Fig. 4, bottom), surface temperature variability is predominantly seasonal (±1.8°C), but cyclic changes associated with ENSO events and a trend component reflecting longer-term warming and cooling are also part of the time series. The partitioning of variability among seasonal, cyclic, and trend components throughout the study area is listed in Table 1. Variability, expressed as standard deviation, is listed as a percentage of the sum of the standard deviations of all components (including error).
For the four pycnocline variables, the state-space decomposition partitions variability of the 5° square time series into components as 35%–56% seasonal, 26%–37% cycles, and 6%–12% trend (Table 1). Stratification differs from the other three variables in having a lower seasonality on average, with a magnitude comparable to interannual variability, and a higher proportion of variability in the trend component.
Trends for the four variables tended to be either monotonic, or with one or two turning points (change of sign of trend) over the 51-yr time series (Figs. 5–8). For some squares, as described in section 2, higher-frequency variability remained in the trend, as can be seen in the gray plots of individual time series trends such as cluster 2 in Fig. 7. This “noise” may have affected the clustering, but we do not feel that it compromised the general patterns of trends presented below.
a. Surface temperature
Surface temperature trends clustered with some regional consistency (Fig. 5). In the eastern equatorial Pacific and eastern tropical Pacific (ETP) warm pool (cluster 1), as well as in parts of the Peru Current and easternmost southern subtropical gyre, surface temperature warmed by ~0.2°C in the late 1970s, but otherwise changed very little. A similar trend occurred at other points in these regions (cluster 2), but a ~0.2°C cooling preceded a warming of ~0.5°C from the mid-1970s to the mid-1990s. In the California Current (cluster 3), along with part of the northern subtropical gyre, surface temperature cooled by ~1°C from 1958 to the early 1970s, then warmed by ~1°C until the mid-1990s, and then cooled again through 2008. In the southern subtropical gyre (cluster 4), surface temperature warmed irregularly by ~0.6°C, mostly after 1980. In the North Pacific (cluster 5), surface temperature cooled by ~0.7°C from the mid-1960s to the mid-1980s and then warmed by ~0.4°C through 2008.
b. Pycnocline depth
Pycnocline depth trends clustered with less regional consistency than for surface temperature (Fig. 6). In the eastern equatorial Pacific (cluster 1), as well as in a few scattered squares in the subtropical gyres and southernmost Peru Current, the pycnocline shoaled monotonically by ~12 m from the early 1960s to 2008. In the ETP warm pool (cluster 2), the pycnocline deepened by ~7 m from the early 1960s to the mid-1980s, then shoaled by ~5 m through 2008. In about half the California Current squares (cluster 3), as well as in some scattered squares in both subtropical gyres, the pycnocline deepened monotonically by ~4 m from 1958 through 2008. In the Peru Current and the eastern North Equatorial Countercurrent (cluster 4), the pycnocline shoaled by ~5 m into the early 1970s, then deepened by ~7 m until the late 1980s, and then shoaled by ~11 m through 2008. In the North Pacific (cluster 5), the pycnocline shoaled erratically until the mid-1980s and then deepened by ~12 m through 2008 with the rate of change accelerating after 2000. In general, the subtropical gyres did not show regionally consistent trends, as indicated by a single cluster corresponding to a coherent region within the study area, and the trend component in many subtropical squares was constant.
c. Mixed layer depth
Mixed layer depth trends clustered with more regional coherence than for pycnocline depth (Fig. 7), in most cases similarly to the pycnocline depth trends. In the eastern equatorial Pacific and North Equatorial Countercurrent (cluster 1), as well as in the southern subtropical gyre, the mixed layer shoaled monotonically by ~6 m from 1958 to 2008; pycnocline depth shoaled by ~12 m in this region. In the ETP warm pool (cluster 3), along with much of the northern subtropical gyre and California Current, the mixed layer deepened monotonically by ~5 m, although most of this change occurred during 1970–2000. A similar deepening trend occurred at other points in these regions (cluster 2), but most of the change in this cluster occurred during 1960–80. The similarity of trends in these two clusters is likely an artifact of our choice to specify five clusters for all the variables in the dataset. The mixed layer depth deepening in these two clusters was similar to the pycnocline deepening in the California Current and subtropical gyres (Fig. 7, cluster 3). In the Peru Current (cluster 4), mixed layer depth was relatively constant or deepened slightly through the 1990s, after which it shoaled by ~10 m; this trend was similar to pycnocline depth in the same region (Fig. 7, cluster 4) in that there was a slight deepening in the 1970s and 1980s and a greater shoaling after 1990. In the North Pacific (cluster 5), the mixed layer shoaled by ~3 m from 1958 to about 1980 and then deepened by ~10 m, again consistent with the pycnocline depth trend in this region.
Stratification trends also clustered with some regional coherence (Fig. 8). In the eastern equatorial Pacific (cluster 1), stratification decreased monotonically by ~50 J m−3 from 1958 to 2008. Both in the ETP warm pool and subtropical gyres (cluster 2) and in the NECC and California Current (cluster 3), stratification decreased by 30–50 J m−3 from 1958 to the early 1970s and then increased by ~80 J m−3 from the mid-1970s to mid-1990s, after which it decreased by 40–80 J m−3. The major difference between trends in clusters 2 and 3 is that the 1970s to 1990s increase in stratification occurred more abruptly in cluster 2. Additional squares in these regions (cluster 4) showed a similar, but much less pronounced, trend of decrease through the 1960s and increase to the mid-1990s, followed by a decrease. In the North Pacific (cluster 5), trends in stratification were more irregular but were marked by an increase of ~100 J m−3 from the late 1980s to late 1990s, followed by a similar decrease to the mid-2000s. This trend had a similar pattern to stratification in the California Current (cluster 3).
Pycnocline depth and stratification trends for 1958–2008 can be summarized by region as follows (stratification changes are expressed as percent of the study area median, ~1000 J m−3):
California Current: Depth varied on decadal scales in some squares but overall deepened by ~5 m from 1970 to 1990s; no overall change in stratification but weakened 5% to the early 1970s, then strengthened 8% to the mid-1990s, and then weakened 4%.
Equatorial Pacific: Pycnocline shoaled by 10 m and weakened 5% overall.
ETP warm pool: Little overall change with decadal-scale variations; deepened 8 m from early 1960s to mid-1980s, then shoaled 5 m; weakened 3% to early 1970s, strengthened 8% during late 1970s, and then no change until the late 1990s when it weakened 8%.
North Pacific: Little change in depth until late 1990s, then deepened by 12 m; little overall change in stratification, weakened 8% during 1960s, then strengthened gradually until abrupt 9% strengthening during 1990s, and then weakened by 10%.
Peru Current and subtropical gyres: No regionally consistent trends, but a tendency to decadal-scale variations with trends reversing every 10 to 30 years.
Table 2 briefly summarizes regional trends for all four variables.
5. Discussion and conclusions
Temporal variability of 1958–2008 monthly sea surface temperature, pycnocline and mixed layer depths, and upper-ocean stratification in the eastern tropical and northern Pacific is primarily interannual in equatorial and tropical waters, and seasonal at higher latitudes. ENSO dominates physical and biological variations in the eastern tropical Pacific (Wang and Fiedler 2006) and its impacts are known to affect the global ocean–atmosphere system (McPhaden et al. 2006). In the state-space model applied to the time series here, the seasonal and cycle components are stationary; these are dynamics that did not change during 1958–2008. Trends represent changes in the dynamics of the variables. For example, Mendelssohn et al. (2005) showed that the ENSO events in several climate indices (Southern Oscillation index, Niño-3 SST, Northern Oscillation index) were captured by a stationary cycle, but that trend terms described changing dynamics that appear as more frequent El Niño events since the 1950s. We observed lower-frequency, long-term trends that showed spatially coherent patterns of variations. Some of these trends were monotonic over the 50 years of this study whereas others showed important decadal-scale variations but little net change. Both forms of long-term variability may be related to global climate change and have ecological implications. The observed surface temperature and pycnocline trends tended to be more regionally consistent in the California Current, North Pacific transition zone, and eastern equatorial Pacific and less so in the eastern Pacific warm pool and subtropical gyres.
Surface temperature trends showed an overall warming in the eastern tropical Pacific and subtropical gyres, but overall cooling in the California Current and North Pacific. Deser et al. (2010) found positive linear SST trends over most of the global ocean for several 1900–2008 datasets, but they showed that interannual and multidecadal variability was also present in these SST time series. Empirical orthogonal functions (EOFs) are an alternative, and computationally simpler, way to analyze spatiotemporal variability. EOF analyses of global SST variability (Enfield and Mestas-Nuñez 2000; Messié and Chavez 2011) show a strong influence of multidecadal modes related to the Pacific decadal oscillation (PDO; Mantua and Hare 2002) and North Pacific Gyre Oscillation (NPGO; Di Lorenzo et al. 2008) in the North Pacific and California Current. These modes resemble the corresponding regional trends in Fig. 6 (e.g., correlations with the PDO are +0.48 for the cluster 3 median trend and −0.49 for the cluster 5 median trend). The California Current cooled from 1958 to the mid-1970s, warmed through the mid-1990s, and then cooled again. The North Pacific warmed from 1958 to the mid-1960s, cooled to the mid-1980s, and then warmed slightly. Thus, the North Pacific and California Current portions of the study area are somewhat out of phase, which is a defining characteristic of the Pacific decadal oscillation. Regional differences in long-term variability may be important for future climate changes. Global warming projections (e.g., Solomon et al. 2007, their Fig. 10.8) show some regional variation in projected surface temperature increase over the next century: 2° to 2.5°C over most of our study area but greater in the North Pacific and equatorial regions and less in the southern subtropical gyre and Peru Current regions.
Pycnocline and mixed layer depth trends were very similar, although pycnocline depth changes were greater than mixed layer depth changes. This is expected since the bottom of the mixed layer as defined here is the top of the pycnocline. Pycnocline and mixed layer depth trends showed an overall shoaling in much of the eastern tropical Pacific but an overall deepening in the California Current. Regionally consistent trends were not found in the North Pacific and subtropical gyres, oceanographic systems that were only partially covered by the study area.
Stratification should be closely related to both surface temperature and pycnocline depth. As surface waters warm, either seasonally or longer term, stratification may increase, although other factors such as precipitation, wind mixing, or upwelling may also affect stratification. In the southern California Current system, Kim and Miller (2007) found that surface temperature cooled from 1950 to 1976, warmed from 1977 to 1998, and then cooled again, and that stratification, but not thermocline depth, changed concomitantly. These surface temperature and stratification trends are similar to those observed here for the California Current (Figs. 5 and 8). Similar to the lack of consistency we observed for some variables in some regions, Palacios et al. (2004), in a state-space analysis of variability of stratification at well-sampled locations in the California Current, found onshore–offshore and north–south differences in patterns of long-term variability. Such variation within a region may be the result of spatial variations in local and remote forcing, as demonstrated for coastal upwelling in this region by Macias et al. (2012).
Long-term or decadal-scale changes in the depth and stratification of the pycnocline will have significant effects on primary production and ecosystem processes. High biological production in coastal upwelling ecosystems such as the California and Peru Current systems, and in the equatorial upwelling system of the eastern Pacific, is driven by the introduction of nutrients from the deep ocean in response to wind forcing. Pycnocline depth and stratification in the water column regulate the entrainment of nutrients from depth, such that an anomalously deep and strong pycnocline, if persistent over long periods, can result in long-term reductions in biological production and in ecosystem changes. Le Quéré et al. (2003) demonstrated how climatic changes in ocean stratification can affect phytoplankton production in both nutrient- and light-limited systems. However, in deeply mixed subtropical gyres, additional factors such as changes in local wind and buoyancy forcing or advective nutrient fluxes may also be important (Lozier et al. 2011).
We observed a deepening of the pycnocline in parts of the California Current since 1970 (Fig. 6, Table 2), along with decadal stratification changes throughout the region (Fig. 8, Table 2). The trend toward a deeper, more stratified pycnocline from the mid-1970s to late 1990s is consistent with long-term fluctuations in biological production within the California Current that have been associated with climatic variability on decadal scales (Chelton et al. 1982; McGowan et al. 2003). McGowan et al. (2003) documented a 17% deepening in the depth of the 12°C isotherm (a proxy for thermocline depth and nutricline depth) for the period 1976–2000 with respect to 1950–75, in the California Current off southern California. Coincident with this deepening was a 74% decline in zooplankton displacement volume, indicating a change in zooplankton community composition (Lavaniegos and Ohman 2007). These authors also reviewed a compelling body of evidence pointing to dramatic declines in marine populations involving at least four trophic levels and a large number of taxa over the same period. The downward trend in zooplankton volume in this region has continued through 2009, although it was interrupted by a step increase following the 1997–98 El Niño event (Fig. 18 in Bjorkstedt et al. 2010).
In the eastern equatorial Pacific, we observed a 10-m shoaling of the pycnocline since the mid-1960s (Fig. 6, Table 2), while stratification decreased overall east of 140°W (Fig. 8, Table 2). These changes would be expected to increase nutrient supply to surface waters and the nutrient-limited production of phytoplankton biomass. A global assessment of recent trends in primary production based on satellite data show a general increase in this region, although the trend is not significant for a large fraction of the region (Henson et al. 2010). North of the equator, in the ETP warm pool, we observed decadal-scale changes in pycnocline depth and stratification. The changes in these two variables were uncorrelated and resulted in little or no net change over 50 years: deepening and weakening from 1960 to early 1970s, then deepening and strengthening to early 1980s, then shoaling and no change in stratification to late 1990s, and then shoaling and weakening to 2008. These decadal-scale variations may have affected primary production and ecosystem processes in the ETP warm pool; the Henson et al. (2010) analysis shows an increase in satellite-derived productivity in this region during 1997–2007, consistent with the shoaling pycnocline. In general, ENSO-scale interannual variability is predominant in this region (Fiedler 2002b).
Climate variability can affect marine organisms by individual-level changes in physiology or behavior, population-level changes in life history (mortality, growth, and reproduction), and ecosystem-level changes in productivity and food web interactions (Pörtner and Peck 2010). The pycnocline changes observed here can have direct effects on phytoplankton productivity, with potential indirect effects on populations and processes at higher trophic levels. There has been a great deal of historical analysis, model prediction, and speculation about future changes in marine ecosystems driven by global warming. In general, warming will either increase or decrease primary productivity depending on the local effects on stratification, kinetics of photosynthesis, light limitation, and a myriad of other processes (Sarmiento et al. 2004; Brander 2010). Effects at higher trophic levels are even more difficult to predict when effects on life history processes are considered (Drinkwater et al. 2010). For example, King et al. (2011) synthesized knowledge on physical and population processes in the California Current ecosystem to produce empirical scenarios of expected changes in representative fish species resulting from climate change. This effort is beyond the scope of the current study. We have shown that significant long-term variations or trends in regional pycnocline characteristics have occurred during the past 50 years in the eastern tropical and North Pacific Ocean. These trends underlie seasonal and interannual (ENSO) variations that are also known to affect phytoplankton production and ecosystem processes, and are likely to have consequences for regional populations and ecosystems.
We thank Sam McClatchie and Franklin Schwing for comments on an earlier draft of this manuscript. Three anonymous reviewers provided valuable comments that greatly improved the paper.