Trends in sea surface temperature off the coast of Ecuador and the major processes that contribute to them

Journal of Marine Systems 164 (2016) 151–164

Contents lists available at ScienceDirect

Journal of Marine Systems journal homepage: www.elsevier.com/locate/jmarsys

Trends in sea surface temperature off the coast of Ecuador and the major processes that contribute to them Laurence C. Breaker a,b,⁎, Hans R. Loor c, Dustin Carroll d a

Moss Landing Marine Laboratories, United States University of Delaware, United States c University of Manta, Ecuador d University of Oregon, United States b

a r t i c l e

i n f o

Article history: Received 18 March 2016 Received in revised form 25 July 2016 Accepted 4 September 2016 Available online 06 September 2016 Keywords: Sea surface temperature Coast of Ecuador Long-term trends Warming from 1900 to ~1940 Ensemble empirical mode decomposition LOWESS smoothing Pacific decadal oscillation ENSO

a b s t r a c t Monthly-averaged sea surface temperatures (SSTs) from seven adjoining sub-regions spanning the coast of Ecuador have been examined for long-term trends and the major processes that contribute to them. The study area extends from the coast out to 84°W longitude, and from 4°S to 2°N latitude. The period of observation extends from 1900 through 2014, a period of 115 years. Linear trend analysis shows that the slopes are positive in all sub-regions and statistically significant in 5 out of 7 of them. Based on the trends SSTs increased from +0.10 °C to +1.42 °C over the entire period with the greatest increases occurring in the sub-regions next to the coast and to a lesser extent in the northern and southern subregions. A second non-linear trend analysis was conducted using Ensemble Empirical Mode Decomposition (EEMD). The results from EEMD for the various sub-regions indicate that SSTs increased significantly from 1900 to at least 1940 in the majority of cases. After 1950, SSTs tended to fluctuate with a minimum that often occurred during the 1960s or 1970s. Since circa 1990, SSTs have experienced little change. Because the data tended to be spatially homogeneous, with two sub-regions being possible exceptions, a mean data set was calculated and likewise subjected to EEMD. The results show that SSTs have increased by approximately 1 °C since 1900 but the process has not been uniform. The increase in temperature between 1900 and at least 1940 approaches 1 °C, similar to the trends found in the majority of sub-regions. Between ~1950 and 1990, temperatures decreased slightly until the mid-1960s and then gradually increased until about 1990. Since ~1990 there has been no significant change in SST. Overall, the results of the global mean analysis are generally consistent with the sub-regional results. Several sources of variability have been tentatively identified that may contribute to the long-term changes in temperature that we have observed. Both spectral analysis and EEMD suggest the importance of the Pacific Decadal Oscillation as possibly the greatest contributor. The two major maxima in the PDO index that occurred during the past century (~1940 and ~1990) are generally consistent with maxima and/or inflection points that occur both in the sub-regional and regional long-term trends. It has been difficult to associate specific ENSO episodes with the long-term trends. However, ENSO and the PDO are closely related. Recent modeling studies have shown that the PDO essentially owes its existence to ENSO. Thus, ENSO appears to be a major contributor to the long-term trends through its primary contribution to the PDO. Finally, spectral analysis also revealed the existence of an oscillation whose period is ~73 years, easily of sufficient length to influence the long-term trends. Climate-related oscillations in the range of 70–100 years are often associated with the Gleissberg sunspot cycle and this could be a contributing factor. The long-term trends are not necessarily a dominant feature in our data. However, had the observations been limited to the period from 1900 to ~ 1940, the trends would have been, without a doubt, dominant. Overall, the likely influence of the PDO on the observed long-term trends in SST must be emphasized. Finally, over the past 25 years there has been no significant change in surface temperature over the study area and this may be due in part to the cooling influence of the PDO since ~1990. © 2016 Elsevier B.V. All rights reserved.

1. Introduction ⁎ Corresponding author. E-mail address: [email protected] (L.C. Breaker).

http://dx.doi.org/10.1016/j.jmarsys.2016.09.002 0924-7963/© 2016 Elsevier B.V. All rights reserved.

In this study we examine long-term trends in Sea Surface Temperature (SST) off the coast of Ecuador and the major processes that

152

L.C. Breaker et al. / Journal of Marine Systems 164 (2016) 151–164

contribute to them. The region of interest extends from the coast of Ecuador out to 84°W, and from 4°S, to 2°N. The study area is divided into seven sub-regions. Each sub-region is 2° of latitude by 2° of longitude. Fig. 1 shows a map of the study area and its sub-regions or boxes. Within each box we have acquired monthly-averaged SSTs from 1900 through 2014. This is an unusually long record and may be the longest record of sea surface temperature to be examined off the coast of Ecuador. We note that although there has been extensive oceanographic data acquired and research conducted in the Equatorial Pacific between the western Pacific (at least as far west as 130E), and the Galapagos Islands since the 1970s, less interest has been shown in the region that lies between the Galapagos Islands and the coast of Ecuador. This observation is consistent with an internet search for oceanographic analyses specifically addressing the region east of the Galapagos where a dearth of relevant information was found. We were motivated, in part, to conduct this study to see if obvious or even not so obvious indications of recent global warming could be found in these data that might give support to other studies of a similar nature that have been conducted in the eastern tropical Pacific. Finally, the study area encompasses one of the most dynamic and productive marine ecosystems in the world (Martinez-Ortiz et al., 2015) where SST is expected to serve as a key indicator of fish behavior. The study area comprises an oceanic region influenced by the Peru or Humboldt Current which transports cool waters across the southwest corner of the study area as it travels from the coast of Peru to the Galapagos Islands. Flow near the coast of Ecuador and south of the Gulf of Guayaquil is weak (b ~10 cm/s) and generally to the northwest, more-or-less paralleling the coast (Wyrtki, 1965). Further offshore the flow intensifies and becomes more westerly as the influence of the Peru Current starts to reveal itself. Off the coast of northern Ecuador, the relatively warm Panama Current flows in a southwesterly direction as it approaches the Equator before turning westward toward the Galapagos. This current enters the study area across its north face. Overall, the mean annual circulation is generally weak but increases near the western boundary of the study area where it becomes westerly at speeds of 5–10 cm/s (Kessler, 2006).

The Equatorial Undercurrent (EUC) shoals in the eastern tropical Pacific and as a result contributes cooler waters to the surface layer as it passes through the Galapagos Islands on its path eastward. According to Lukas (1986), the influence of the EUC extends to the coast of Ecuador. These waters provide a cooling influence on SSTs in the study area near the Equator. During ENSO episodes, however, the EUC is driven downward by warmer waters from the western Pacific and its influence at or near the surface is greatly reduced. The processes that affect SST and their interaction in the tropics form a complicated picture. According to Sarachik (1978), equilibrium surface temperatures near the Equator are determined by the interaction between short-wave solar radiation, infrared radiation, and evaporation, and so are strongly influenced by atmospheric processes. More specifically, with regard to the eastern tropical Pacific, SST tends to be a balance between four processes (Behringer et al., 1998). First, is the strength of upwelling along the Equator which is determined by the westerly component of the wind. Stronger winds to the west produce greater Ekman transport away from the Equator leading to increased upwelling. Second, is the northwesterly transport of colder waters from the coasts of Peru and Ecuador toward the Galapagos. Third, is the mixing of cooler Equatorial waters with warmer waters to the North and South of the Equator. Finally, the fourth process is heat exchange between the ocean and the atmosphere along the Equator. When changes from non-El Niño to El Niño conditions occur each of these processes is strongly affected. The study area has two rather distinct oceanic seasons: summer, which occurs between January and April, and winter, which occurs between July and October. Other times of the year are considered to be transition periods between summer and winter (Cucalon, 1989). The study area is also a transition zone that separates two water masses, Tropical Surface Water to the north, and Subtropical Surface Water to the south. Between these water masses lies a major oceanic boundary called the Equatorial Pacific Front (Laevastu and LaFond, 1970). This frontal region is characterized by strong property gradients including SST. According to Cucalon (1989) this front is unpredictable during the summer regarding both its location and strength. During the winter

Fig. 1. This figure shows the oceanic region employed in this study. Monthly-averaged Sea Surface Temperatures (SSTs) from 1900 through 2014 provide the basis for this study. The region is subdivided into seven sub-regions or boxes and are numbered accordingly. Each box covers an area that is 2° of latitude by 2° of longitude (2° × 2°). Boxes closest to the coast (3, 5, and 7) are partially covered by land. The center locations for each box are as follows: Box 1: 1°N, 83°W; Box 2: 1°N, 81°W; Box 3: 1°N, 79°W; Box 4: 1°S, 83°W; Box 5: 1°S, 81°W; Box 6: 3°S, 83°W; Box 7: 3°S, 81°W.

L.C. Breaker et al. / Journal of Marine Systems 164 (2016) 151–164

the front becomes more well-established and thus easier to identify. Equatorial upwelling has been reported east of the Galapagos extending to the coast of Ecuador, and may be associated with the Equatorial Pacific Front in this region, according to Pak and Zaneveld (1974). Equatorial Kelvin waves that are generated in the western Pacific and propagate eastward enter the study area as a result of El Niño. These waves reflect off the coastal boundary and propagate poleward along the coasts of North and South America as coastal Kelvin waves (Clarke, 1983). Associated with equatorial Kelvin waves are eastward flowing currents, often called El Niño currents that transport warmer waters toward the coast of Ecuador often increasing water temperatures by several degrees Centigrade. During the 1982–83 El Niño, SSTs along the coast of Ecuador rose by as much as 6.3 °C (Cucalon, 1987). These episodes occur approximately every 3 to 7 years, often lasting for up to a year (e.g., Philander, 1990). Major El Niño events, which occur every 10–15 years, can last for more than a year. The envelope formed by connecting the major events during the past and present centuries is referred to as the El Niño – Southern Oscillation (ENSO) intensity modulation cycle (e.g., Sun and Yu, 2009). Interaction between the Pacific climate system and ENSO may be responsible for the modulation cycle. A study published earlier by Breaker et al. (2015) provides a preliminary analysis of the data used in this study. They found that SSTs generally increase toward the coast and decrease toward the south over the same region employed in this study. However, the variability in SST decreases toward the coast and increases toward the south. El Niño episodes, particularly those that occurred in 1982–83 and in 1997–89, were dominant features in each record. Trend analyses showed that surface temperatures, overall, increased between 1900 and 2015 at each location and that in most sub-regions major increases in SST occurred between 1900 and at least 1940 at rates approaching 1 °C/50 years. Between about 1950 and 1990, SSTs decreased slightly and then increased. After ~1990, surface temperatures varied only slightly with no obvious trends during this period. The first study was preliminary in nature and primarily descriptive. As part of this earlier work the long-term trends were briefly examined. In this study we begin with a more detailed analysis of the long-terms trends followed by an examination of the processes that are most likely responsible for producing them.

2. The data As indicated in the introduction, the region of this study extends from the coast of Ecuador out to 84°W, and from 4°S, to 2°N. The study area is divided into seven sub-regions or boxes that are each 2° × 2° of latitude and longitude (Fig. 1). The spatial distribution of the sub-regions is shown in Fig. 1. Within each sub-region monthlyaveraged SSTs have been acquired for the period from 1900 through 2014, yielding a record that is 115 years in length. These data are from the International Comprehensive Ocean-Atmosphere Data Set (ICOADS) and were provided by NOAA/OAR/ESRL PSD, in Boulder, Colorado, USA.

153

In more detail, the individual observations within each box, regardless of their location within that box for a given month, have been averaged to provide the monthly mean values. Locations for the original observations or the exact times they were collected are not available. Thus, we cannot assign exact locations for the monthly values within any box. This restriction plus the fact that there are only seven boxes limits the possibilities for conducting meaningful spatial analyses such as Empirical Orthogonal Function analysis, and is the reason why we have focused on the time domain. The number of observations available within each sub-region is given in Table 1. Further details concerning the quality-control procedures that were employed and the data processing techniques that were used can be found at the following website: http://www.esrl.noaa.gov/psd/. After receiving the data from ICOADS, they were initially screened to extract the monthly mean values of SST for each box and the number of observations that entered into the calculation of the mean values for that month, together with the corresponding months and years. These values were written to seven separate files, one for each sub-region that were subsequently used to conduct the various analyses. At this point we applied a correction to the data based on the work of Folland and Parker (1995). This correction applies to data between 1856 and 1941 and is related to the use of canvas and wooden buckets to store the sea water samples from which the measurements of SST were subsequently obtained. During this period a positive bias correction was derived that starts at 0.11 °C in January 1856 and ends at 0.42 °C in November 1941. This increasing bias was due to the increasing use of canvas buckets and increasing ship speeds over the entire period. We do question the abrupt cutoff for this correction in November 1941 since it assumes that all ships started using the correct procedures at the same time when in reality it probably took several years for these improvements to take effect. In any case, the correction is approximately linear and we have assumed a linear relationship between the correction and the period in question to obtain the corrections we have applied to the data. The results are shown in Fig. 2. The upper panel (Fig. 2a) shows a comparison of the corrected and uncorrected data for a brief period just prior to the cutoff date after which the correction was no longer applicable. The correction increases the values of SST by a small but noticeable amount. In Fig. 2b, we include the entire record and show the data together with smoothed versions of same. A close look at the smoothed versions of the data show the gradual increase in the correction from 1900 to about 1930, beyond which the termination of the correction and the smoothing that was applied reduces the difference between the two curves as they approach the same value. At this point, it becomes apparent that the correction is significant and its absence could affect the interpretation of our results later in the study. Data density is an important issue since the number of observations per month are limited during the early years of this study and during WWII. Fewer observations during these periods translates directly into greater uncertainty in any analysis we conduct with them. Plots of the number of observations per month for each sub-region are shown in Fig. 3. To see more clearly what the data densities were particularly prior to say 1950, we have plotted only the number of observations

Table 1 Selected properties of the data and data analysis parameters for each sub-region. The underlined values in row 4 indicate that the slopes are not statistically significant at the 95% level of significance.

Number of Observations (% of 1380 months) Sample mean (°C) Sample variance (°C2) Linear Slope (°C/year) Overall (linear) change (°C) Values of Nstd

1°N, 83°W

1°N, 81°W

1°N, 79°W

1°S, 83°W

1°S, 81°W

3°S, 83°W

3°S, 81°W

Regional mean values

988 (71.6%) 25.71 1.80 +0.0076 +0.87 0.48

1269 (91.9%) 25.97 0.90 +0.0088 +1.01 0.45

953 (69.1%) 26.49 1.25 +0.0123 +1.42 0.42

877 (63.6%) 24.14 3.19 +0.0009 +0.101 0.38

1264 (91.6%) 24.81 1.75 +0.0022 +0.25 0.38

859 (62.2%) 22.28 7.43 +0.00944 +1.08 0.34

1286 (93.1%) 23.38 3.29 +0.0129 +1.42 0.41

N/A 24.68 1.28 +0.0059 +0.674 0.28

154

L.C. Breaker et al. / Journal of Marine Systems 164 (2016) 151–164

Fig. 2. Annual global average positive bias correction applied to the regional mean record of sea surface temperature for the period between 1900 and 1941. Fig. 2a shows a detailed view for the last 3 years of the correction period with the original data in black and the data with correction in green. Fig. 2b shows the impact of this correction over the entire period. In this case the red (corrected) and blue (original) curves have been smoothed using LOWESS with α = 0.11. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

per month up to a maximum of 100 (values N ~100 are obviously of less concern and are not shown). This permits us to look carefully at the data available during the early years. In a few cases during the first decade no observations were available for certain months. The number of observations per month in each 2° × 2° box since 1900 was often b25 until at least the 1960s. Although it is not clear why, we note that for the three sub-regions at 81°W there is a significant increase in the number of monthly observations between about 1922 and 1932. Major gaps occur in each record in the mid-1940s during WWII when routine data collection at sea was essentially discontinued. From 1944 through 1946 the data are sparse and typically many months are missing during this three-year period. By the 1970s, the numbers of observations in

each sub-region have generally increased to 100 or more per month almost certainly due in part to the availability of satellite data. Although boxes 3 (1°N, 79°W), 5 (1°S, 81°W), and 7 (3°S, 81°W) are partially covered by land, Table 1 shows that they still contain more observations than several boxes that are land-free (boxes 4 and 6). Because of the gaps, the data are non-uniformly spaced in time and so interpolation was often required to produce time series that were uniformly spaced. For cases where uniform spacing was required, such as in EEMD, we have used a nearest-neighbor approach in the earliest years where we have repeated the nearest 12 monthly values to construct a continuous record that reproduces the annual cycle back to 1900. This procedure was required in a few cases over a period of only

Fig. 3. This figure shows the number of observations per month for each sub-region for values of less than or equal to 100. This provides a more detailed examination of the data density during the early years of the record and WWII where the data densities are sparse.

L.C. Breaker et al. / Journal of Marine Systems 164 (2016) 151–164

one or two years between 1900 and 1902. For missing data further into the records where no extrapolation was required we have used a cubic spline to fill the gaps (e.g., Davis, 1986). In any case, Table 1 shows that even in the worst case (box 6), 62% of the original data exist. However, as we continue to examine these data we have gained confidence in our results due to the fact that the trends we observe are well-defined and generally consistent between the different sub-regions. What is fortunate in many cases is that where the data were missing the gaps were often relatively short. Finally, the issue of data density has influenced not only the methods of analysis we have chosen to employ but the goals we have chosen for this study as well. We are looking primarily at long-term trends where gaps are far less of a problem than it would be if we looking a higher frequency variability. It has been our experience that the data can be seriously undersampled and the underlying trends and related lowfrequency variability of primary interest can still be extracted without serious problems. To illustrate this point, we return to Fig. 2b, where the smoothed corrected and uncorrected curves illustrate not only the differences we discussed due to the correction that was applied but also they illustrate the ability of our data to resolve all of the major and some of the lesser El Niño episodes since 1900 (1917, 1925–26, 1932, 1940–41, 1957–58, 1972–73, 1982–83 and 1997–98). 3. The approach In this study we employ several methods of analysis to (1) extract trends in the data, and (2), identify oceanic processes that contribute to these trends. Three methods of trend analysis are used, linear leastsquares trend estimation, Ensemble Empirical Mode Decomposition (EEMD) and the LOWESS smoother. The first method is well-known and so will not be elaborated upon (e.g., Natrella, 1963). The second method, EEMD, is also well-known and has been expounded upon frequently (e.g., Huang et al., 1998; Huang, 2005a,2005b). This method decomposes a time series into a sequence of empirically orthogonal Intrinsic Mode Function components (IMFs) and a residual (Huang et al., 1998). In EEMD, the number of modes is determined by the data. Thus, EEMD is data adaptive, and is also well-suited for the analysis of non-stationary and nonlinear time series. The IMF components are often physically meaningful because the characteristic scales in each case are determined by the data. However, selected modes may require grouping in order to extract a physical basis. In practice, extracting the trend, depending on the definition that is employed, may require that several of the higher modes including the highest mode (the so-called residual) be grouped to reconstruct it. This selection or grouping, off course, is strongly influenced by the time scales of variability that we wish to include in the long-term trends. Finally, the method of extracting the individual modes has been improved by including a noise-assisted component in its calculation (Wu and Huang, 2009), which defines the true IMF as the mean of an ensemble of IMFs, and, in the process, preserves the physical uniqueness of the decomposition. As a result, the name of the method was changed from simply Empirical Mode Decomposition (EMD) to Ensemble Empirical Mode Decomposition (EEMD). Because of this modification two parameters must now be specified, a noise-related parameter that is obtained from the data initially, Nstd, and the size of the ensemble that is employed. Nstd is calculated as the ratio of the standard deviation of the first mode in the decomposition to the standard deviation of the original data. The values of Nstd used for each sub-region are included in Table 1. In this study we have used an ensemble size of 300 in all cases. The third method is called LOWESS and is a state-of-the-art smoothing operator that can be used to extract meaningful patterns that may not be apparent in the original data (Cleveland, 1979, 1985). LOWESS is a non-parametric technique for performing robust, locally-weighted regression. The method fits linear or quadratic basis functions of the predictor at the center of neighborhoods where the radius of each

155

neighborhood contains a specific percentage of the data points. This type of smoothing function introduces minimal distortion at the boundaries of the data. The fraction of data in each neighborhood, and thus the smoothness, is determined by a smoothing parameter specified by the user. In the application of LOWESS, two parameters must be specified, the degree of the local polynomial that is employed (linear or quadratic) and the level of smoothing, α. Throughout this work the degree of the local polynomial has been set to quadratic. The smoothing parameter is specified by the choice of α, where 0 ≤ α ≤ 1. For α equal to 0, no smoothing occurs and for α equal to 1, we obtain essentially a straight line. The goal is to choose α large enough to obtain as much smoothness as possible without distorting underlying patterns in the data. To assist in identifying oceanic processes that contribute to the trends in SST, we employ two methods of spectral analysis, the Maximum Taper Method (MTM) of power spectrum analysis (Thompson, 1982) and the periodogram (e.g., Chatfield, 1995). The MTM method is widely used and is considered to be state-of-the-art in spectral analysis. It provides an improved estimate of the power spectral density by employing a bank of optimal bandpass filters to provide the smoothing. It also provides a variable time-bandwidth product that yields a balance between the variance and resolution. The periodogram is a timehonored way of achieving high spectral resolution at the expense of lower confidence in the spectral estimates. Together, the MTM method and the periodogram provide and excellent way of achieving a balance between confidence in the spectral estimate and the resolution. The final technique we employ is cross-correlation analysis. In the field of time series analysis, cross-correlation provides a measure of the similarity between two times series as a function of the lag between one series and the other. A complete discussion of this function and its calculation can be found in Jenkins and Watts (1968). In addition, we have calculated the cross-correlation functions with pre-whitening. Pre-whitening removes serial or autocorrelation from a time series and is a recommended procedure in conducting cross-correlation analysis of climatological data (Katz, 1988). Pre-whitened cross-correlations are invariably lower in magnitude than their non-pre-whitened counterparts because the confounding effects of autocorrelation within each series have been at least partially removed. Further details can be found in Weedon (2003). 4. The results 4.1. Preliminary results First, we examine the monthly SST data shown in Fig. 4 for each subregion. The data have been lightly smoothed (yellow) using LOWESS (α = 0.045). Even at the scale shown we can observe certain features of interest. Major El Niño episodes that occurred in 1982–83 and in 1997–98 can most readily be observed at 1°S and 3°S suggesting an increased response to El Niño warming events moving southward along the coast of Ecuador. We also observe that the variability in SST between the 7 sub-regions varies significantly. This variability tends to increase in sub-regions centered at 1°S, 83°W, 3°S, 83°W and 3°S, 81°W. This increase in variability is related in part to the increasing influence of the Southeast Trade Winds further to the south and offshore which intensify between May and December. We can also see a tendency for increased variability in SST on interannual time scales within each subregion that appears to be associated with El Niño events. During these episodes the equatorial wind regime reverses in direction from essentially eastward to westward. Finally, the sparsity of data during WWII (~1944–1946) can readily be observed in several sub-regions including 1°N, 79°W, 1°S, 81°W and 3°S, 81°W. In this section we also refer the reader to Table 1. The number of observations of SST for each sub-region are shown in row 2 of the table. The mean values for each sub-region are shown in the third row of the table and the variances for each sub-region in the fourth row. The linear slopes for each sub-region are shown in row 5, the overall change

156

L.C. Breaker et al. / Journal of Marine Systems 164 (2016) 151–164

Fig. 4. Monthly-averaged sea surface temperature (°C) for each sub-region. The records start January 1900 and end December 2014. The original data are shown in black. The data have been smoothed using LOWESS with the smoothing parameter set equal to 0.045 (yellow). The vertical scale has been kept constant to facilitate comparisons. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

in temperature from the beginning to the end of each record is shown in the 6th row of the table, and the values used for the noise-related parameter, Nstd, in the EEMD decompositions are shown in row 7. 4.2. Analysis of trends Next we examine the long-term trends in SST from each sub-region. First, we calculate the linear trends. They are all positive but the two sub-regions centered at 1°S, 81°W and 1°S, 83°W are not statistically significant at the 95% level of significance. Correctly interpreted this means that these slopes are equally important but they are not significantly different from zero. We have also tabulated the overall increases in SST at each location based on the differences between the first and last values of the linear fits to the data. (The overall increase or range in SST could have been determined in several ways - we chose to extract the range from the linear approximations). The largest increases in temperature occurred at 3°S, 81°W (+1.42 °C) and at 1°N, 79°W (+1.42 °C). Two patterns emerge: First, the greatest increases in SST occur in the sub-regions next to the coast, and second, greater increases occur in the northern-most and southern-most sub-regions – or, restated, the minimum increases occur at 1°S, 83°W and 1°S, 81°W. Although a linear model was first employed, inspection of the monthly data suggested that linear models might be a poor choice since significant long-term changes occur in all of the records that would not be well represented by a linear fit to the data. As a result, we next employ nonlinear methods in order to identify and extract the long-term trends, and to provide an improved basis for identifying the process or processes that may be contributing to these trends. To allow for non-linear behavior in the data we have applied Ensemble Empirical Mode Decomposition (EEMD) to each record to extract information regarding the long-term trends. The final or last mode in an EEMD decomposition is called the residual and often captures trendlike behavior. In the application of EEMD, the data must be uniformly spaced in time. Because of the gaps that occur we have interpolated

the data to a uniform grid with monthly spacing using the methods described earlier. The lower modes in EEMD contain the higher frequencies and these modes will be influenced by the interpolation, particularly for the longer gaps. As the mode number increases and the corresponding frequencies are lower the influence of the gaps gradually disappears. We estimate that for mode numbers of six and higher (where the oscillations have periods of at least 10 years) the gaps no longer significantly affect the modal patterns. An EEMD decomposition is shown in Fig. 5 for the data from the subregion centered at 1oS, 81oW. This decomposition and the other decompositions each produced 10 modes, consistent with the notion that the data used in this study come from a region that tends to be somewhat homogeneous. The values of the noise parameter, Nstd, that was used in this decomposition and the other decompositions ranged from about 0.28 to 0.48 (Table 1). As stated earlier, the ensemble size was 300 in each case. We note that relatively small changes in either of these parameters usually have only a minor effect on the modal patterns that are produced. The lower modes because of the shorter time scales associated with them are not of interest to us in this study. For mode 8, where nonuniform sampling is not an issue, the oscillatory pattern shown in Fig. 5 appears in the decompositions at the other locations as well. These patterns are characterized by major peaks located around 1940 and 1990, plus or minus a few years. These maxima may correspond to the Pacific Decadal Oscillation (Mantua et al., 1997) which has had a pattern of variability during the last century that is similar to the patterns we have just described. The Pacific Decadal Oscillation (PDO) is a pattern of climate variability that in some respects is similar to ENSO, but has time scales of at least several decades. The PDO, because of its long period and relatively high amplitude, may well contribute to long-term trends in the data. We have also calculated the variance by mode for each sub-region. The results are shown in Fig. 6. The variance or power of the modes tends to decrease with increasing mode number beyond the second mode except for the last mode (i.e., mode 10) where the variances

L.C. Breaker et al. / Journal of Marine Systems 164 (2016) 151–164

157

Fig. 5. An Ensemble Empirical Mode Decomposition (EEMD) of sea surface temperature for the sub-region centered at 1°S, 81°W. The decomposition produced 10 Intrinsic Mode Functions (imfs) or modes. Modes six and higher are considered to be valid since the time scales associated with these modes are long enough to be unaffected by gaps that occur in the data.

Fig. 6. 10Log10 of the variance (i.e., deciBels) plotted versus mode number for each mode from Ensemble Empirical Mode Decomposition for each sub-region. For modes 7–9, the variances in some cases overlapped and so one value was plotted just to the right of the next closest value for clarity.

associated with the residual or long-term trend increase significantly. Also, the variance associated with the sub-region centered at. 3°S, 83°W has the highest value for all of the modes except for the residual. This increase in variability is almost certainly related to greater contact with, and influence from the Humboldt Current where greater variability may be expected. Finally, regarding mode 8, which may be related to the PDO, the variances were significantly higher at 1°S, 83°W and at 3°S, 83°W, suggesting that the influence of the PDO may be considerably greater as we move further offshore. Although our goal is to extract the long-term trend for each subregion definitions for the trend are rare. The term “trend” only acquires meaning when a specific procedure is used in its estimation, according to Fuller (1996). Because subjectivity enters the picture in deciding on what is trend-like and what is not, we have chosen EEMD as a good tool for addressing this problem because the individual modes are determined by the data per se making the process of decomposition to some extent, more objective. However, the problem now does become somewhat subjective because we must decide which, if any modes, in addition the highest mode or residual, should be included as contributing to the long-term trend. This is called grouping and is a crucial step in certain methods of trend extraction where methods of modal decomposition are employed (e.g., Singular Spectrum Analysis). In our case we

158

L.C. Breaker et al. / Journal of Marine Systems 164 (2016) 151–164

wanted the trends to reveal more information regarding the processes that contribute to them than we could obtain from a linear analysis or from the EEMD residual alone. Consequently, we have grouped the highest three modes (modes 8 + 9 + 10) from the EEMD decomposition for each sub-region. Thus, we include time scales as short as interdecadal as part of the long-term trend. As a way to evaluate the EEMD-derived trends, we have also employed the LOWESS smoother. The results from LOWESS provide a basis for comparison and thus evaluation. After some experimentation we adopted a value of 0.45 for the smoothing parameter, α. The results are shown in Fig. 7. The trends from EEMD are shown in red and those from LOWESS in blue. Although the temperatures vary from one subregion to the next, a range of 3 °C has been maintained throughout. Differences between the methods occur during the early years in the subregions at 1°N, 79°W and 3°S, 81°W. Lesser differences can be found elsewhere. We note the differences, although slight, at the ends of the records at 1°N, 83°W and 3°S, 81°W. These differences occur where the trends might provide useful information about the future. To their credit both methods are known to have minimal end effect problems consistent with our own experience in using them. In an effort to determine whether these differences could be related to the methods or to the data, we reversed the records in time and subjected the timereversed data to the same operations. The results were the same in each case. These differences most likely simply reflect slight differences in the frequency response of each method. Had we varied the degree of smoothing in LOWESS or the noise parameter in EEMD we might have been able to obtain even closer agreement in some cases. In the majority of cases, these figures indicate that SST tends to increase significantly from 1900 to 1940 or slightly beyond. After 1950 there is often a well- defined decrease in temperature that occurs between about 1955 and 1980. After ~ 1990 there are no major changes in temperature. In some cases temperature increases slightly during this period and in some cases it decreases. As we will see, when the data from all seven sub-regions are averaged, the slope from approximately 1990 to 2015 is close to zero. All of the EEMD-derived trends are plotted together in Fig. 8, making it easier to observe what these trends have in common. Although these trends do have a number of features in common we note that the trend for the sub-region centered at 3°S, 81°W differs from the rest. The linear

Fig. 8. The long-term trends in SST from each sub-region obtained from EEMD plotted together on the same axes to facilitate comparison. The different colors correspond to the different sub-regions as indicated in the legend. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

trend in this case is slightly positive (Table 1) but the nonlinear trend shown here does not appear to be positive, in contradistinction to the other trends. This location is uniquely influenced by warmer waters which enter the offshore domain from the Gulf of Quayaquil which constitutes a major portion of this sub-region. We return to this issue in the final section of the text. Continuing, the increase in temperature from 1900 to 1940 is very clear in the majority of cases. However, at three locations along the coast (3°S, 81°W, 1°S, 81°W and 1°N, 79°W) we see first a decrease in temperature for at least the first decade followed by increasing or slightly increasing SSTs for varying periods depending on the location. We discuss this behavior in more detail in the final section. A minimum in temperature between 1960 and 1980 can be found in all cases except at 3°S, 81°W where it is more difficult to identify a single minimum. Between the early 1960s and 2000, SSTs increase significantly at every location (~0.3 °C to ~0.6 °C) although the start and end times of increasing temperature vary from one location to another. Finally, at two locations (1°N, 83°W and 3°S, 83°W) increasing SST is observed during the last decade and at 3°S, 81°W decreasing SSTs.

Fig. 7. Long-term trends in SST from Ensemble Empirical Mode Decomposition (red) and the LOWESS smoothing function (blue) plotted separately for each sub-region. The trends from EEMD were obtained by summing modes 8, 9, and 10 in each case. The LOWESS smoothing parameter, α, was set equal to 0.45. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

L.C. Breaker et al. / Journal of Marine Systems 164 (2016) 151–164

are observed during this period. Because these short-term trends are either not consistent with neighboring sub-regions or are of opposite sign (3°S, 83°W and 3°S, 81°W) it is difficult to determine whether or not they are significant. To summarize the trend analyses, the long-term trends in SST from most of the seven regions are roughly in agreement. They can be subdivided into three periods: the period between 1900 and at least 1940 where SSTs generally increase, the period between approximately 1950 and 1990 where SSTs tend to fluctuate, first decreasing and then increasing, and, finally, the period between ~1990 and 2015 where temperatures display no significant trends. 5. Spectral analysis In this and the following section we examine sources of variability in the data that may contribute to the long-term trends. First, we calculated the periodogram for each sub-region and found that the spectra were highly similar (not shown). In each case, the annual cycle at 1 cycle per year, and usually, the first and second harmonics at 2 and 3 cycles per year were the dominant signals in the spectra. Because of the strong similarity in the spectra we have taken an average of the 7 subregions to obtain a mean record that we then subjected to spectral analysis. (We also calculated the mean spectrum from the 7 individual spectra and obtained almost identical results, as we should.) Continuing, we proceeded to calculate the periodogram (black) and the MTM power spectral density (blue) from the mean time series obtained by averaging all 7 records. For the periodogram a rectangular window with a length of 2048 was used to obtain the spectral estimates. For the MTM power spectral densities, a time-bandwidth product of 2 was used to specify the data window whose length was again set at 2048. The timebandwidth product is a variable that can be increased over a small range from approximately 2 to 7 to obtain smoother spectral estimates. Thus, in our case we used less smoothing to obtain higher spectral resolution. The results are shown in Fig. 9. The spectra on plotted on a scale from 0.0 to 1.11 cycles per year (cpy). The vertical axis is logarithmic and plotted in deciBels (dB). The annual cycle at 1 cpy is the strongest signal in the data. We will return to the annual cycle shortly. One other feature in the spectra is that there is generally a negative slope from 0.0 to at least 1.0 cpy. This slope is due to the fact that stronger ocean signals tend to occur at lower frequencies plus the effect of serial correlation or autocorrelation in the data that translates into so-called red noise in the frequency domain which also contributes to the negative slope (e.g., Weedon, 2003).

159

In addition to the annual cycle, frequencies between zero and approximately 0.5 cpy contain the information that is most relevant to this study. As indicated in Fig. 9, in the region between about 0.10 cpy to 0.30 cpy or slightly higher (periods ranging from about 3 to 10 years) we find increased spectral amplitudes that are usually associated with ENSO variability, in this case, in the eastern tropical Pacific. Historically, these events tend to occur every 3–7 years (Philander, 1990). That the higher amplitudes observed in this region (particularly apparent in the blue trace) survived the averaging process further attests to their significance since ENSO activity is not periodic but rather, aperiodic in nature. Next, a significant peak in the spectrum occurs at approximately 0.075 cpy (13–14 years) and may correspond to the ENSO modulation cycle which has a preferred period in the range of 10–15 years (e.g., Sun and Yu, 2009). This modulation cycle connects ENSO maxima of greater intensity that occur less frequently than every 3–7 years. For example, the period between the 1982–83 and 1997–98 M-El Niño events is approximately 15 years. Finally, a major peak in the spectrum occurs at roughly 0.02 cpy with a corresponding period of almost 60 years. This peak may be related to the Pacific Decadal Oscillation (PDO) whose major maxima during the past century occurred circa 1940 and 1990. Although the lack of spectral resolution prevents us from making a positive identification in this case, we are not aware of any other long-period oscillation that has a similar interdecadal time scale. We next examine the mean periodogram in Fig. 9 for side band frequencies that are associated with the annual cycle. Side bands indicate the presence of long-period oscillations that alter the amplitude of the annual cycle in accordance with the modulation period. The influence of modulation is indicated by the presence of well-defined upper and lower side bands that are located symmetrically about the center frequency. Although distinct upper and lower side bands on each side of the annual cycle are apparent in the periodogram they are not resolved in the MTM spectrum. The modulation frequency is given by the difference (in an absolute sense) between the upper or the lower side band frequency, and the center frequency. In this case we estimate the difference between 1 cpy and the associated upper side band frequency because it is better resolved than the corresponding lower side band. The resulting period in this case is approximately 73.0 years. By examining the side bands associated with the annual cycle we have been able to detect an oscillation whose period exceeds half the record length, one that would not have been detectable in the primary spectrum. Such an oscillation might well contribute to the long-term trends we examined earlier. Finally, we examined mode 9 from the previous EEMD analysis (Fig. 5), whose range of periods approach centennial, to see if a similar

Fig. 9. The periodogram (black) and the MTM power spectral density for the mean record obtained by averaging the data from all seven sub-regions. The range of frequencies shown extends from approximately 0.0 to 1.1 cycles per year, and the vertical axis represents the spectral amplitude in deciBels (dB). See text for further details.

160

L.C. Breaker et al. / Journal of Marine Systems 164 (2016) 151–164

oscillation could be identified. The oscillations in mode 9 vary over a rather wide range of periods, and the phases vary significantly from one sub-region to the next making it difficult to find any obvious connection between them. 6. The Pacific decadal oscillation The amplitude, spatial structure and time dependence of the Pacific Decadal Oscillation are related to the Pacific Decadal Oscillation Index which is defined as the leading principal component of monthly SST variability in the North Pacific, poleward of 20°N. However, the influence of the PDO can be found in both the North Pacific and South Pacific Oceans, according to Shakun and Shaman (2009). They have calculated a similar index for the Southern Hemisphere based on SST and found that it adequately explains the spatial structure of the PDO in the North Pacific as well. They conclude that pacific decadal variability (PDV) is similar on both sides of the Equator, suggesting the tropics as a common source for PDV. Thus, although the PDO index is not formally defined for latitudes b 20°N, we expect to find evidence for the existence of the PDO in equatorial regions as well. Returning to Fig. 5, the two maxima that occur in the EEMD pattern associated with mode 8 (imf8), the first around 1940, and the second between 1990 and 2000, most likely reflect the influence of the PDO. We have also taken mode 8 from the other decompositions and plotted them together in Fig. 10. The PDO index has been superimposed (solid red) to show the similarity between the individual modes and the PDO itself. The agreement is not perfect and in most cases there is a lag of several years between the primary peaks in the modal patterns, and the PDO index, consistent with a possible cause and effect relationship. In any case, when the lags are taken into account, the correlations between these patterns and the PDO are highly significant. Finally, we note that since the late 1980s the PDO index has become strongly negative and thus its impact on surface temperature, if it is significant, should be to lower SSTs. Returning to Fig. 8 with continued emphasis on the PDO, there are certain features that the majority of trends have in common. During the period between 1900 and about 1940, the trends are generally in phase with the Pacific Decadal Oscillation index and thus the PDO may have contributed constructively to the positive slopes that occur in most cases during this period. All of the trends show a prominent or less prominent maximum in temperature around 1940 which coincides closely to the first major maximum in Pacific Decadal Oscillation index (Fig. 10). This maximum in SST is usually followed by a minimum that occurs during the 1960s or 1970s, which coincides at least

approximately with the major minimum in the PDO index that occurs between 1955 and 1980. Less apparent is a slight maximum or inflection point that occurs in the temperature trends around 1990 – this feature coincides with the second major maximum in the PDO index. Overall, the apparent influence of the PDO is unmistakable when the PDO index is compared with the SST patterns shown in Figs. 7 and 8. Next, we draw once again upon our results from EEMD to illustrate the likely influence of the PDO in a different manner. In this case we have essentially removed the possible influence of the PDO by subtracting mode 8 from the original data at each location. We can view this operation as that of a band-rejection filter whose center frequency has a mean period for mode 8 of roughly 50 years. However, this filter is far from ideal and so information from adjacent bands is not completely removed. This operation can be expressed as r ij ðt Þ ¼ oij ðt Þ−mij ðt Þ;

ð1Þ

where oij (t) represents the original data, mij (t) represents the data for mode 8, and rij (t) represents the original data less the influence of mode 8, as a function of location, mode number and time, respectively. Specifically, i represents the sub- region, i = 1,2, …, 7, j represents the mode number and is equal to 8 in each case, and t is the time in months from 1 to 1380 months (115 years). The likely influence of the PDO was detected in each sub-region but was found to be somewhat greater offshore than closer to the coast, generally consistent with the variances shown for mode 8 in Fig. 6. We have chosen to show only the results for the sub-regions centered at 1°N, 83°W and 1°N, 81°W (Fig. 11a), and at 1°S, 83°W and 3°S, 83°W (Fig. 11b) where the likely influence of the PDO is most apparent. These figures show the original data (solid curves) and the data with mode 8 removed (dashed curves), smoothed with LOWESS. The smoothing parameter, α, was set equal to 0.42 in each case except for the lower trace in Fig. 11b (3°S, 83°W) where it was set equal to 0.56. In Fig. 11a, in both sub-regions we see that between about 1920 and 1950 the traces with mode 8 removed (dashed) are reduced in amplitude compared to the original data, between about 1950 and 1980 they are slightly higher in amplitude, and between about 1980 and 2000, they are again slightly lower in amplitude than the original data. These differences are consistent with the influence of a cyclic pattern that is similar to the PDO. In the lower panel in sub-regions 1°S, 83°W and 3°S, 83°W (Fig. 11b), we again see that the differences between about 1920 and 1950, and between about 1950 and 1980 are consistent with the influence from the PDO or a PDO-like oscillation whose period is interdecadal. However, between about 1980 and 2000, the differences between the original data (solid) and the data with mode 8 removed (dashed) are slightly smaller than they were in the upper panel (1°N, 83°W and 1°N, 81°W) and although consistent with the influence from a PDO-like oscillation, these differences may not be significant. Although the results obtained be removing mode 8 from the original data are not overwhelming, they are consistent with the notion that if mode 8 does, in fact, portray the Pacific Decadal Oscillation the comparisons above are generally in agreement with our expectations. 7. From sub-regional to regional

Fig. 10. The eighth IMF (imf8) or mode 8 from each Ensemble Empirical Mode Decomposition (dashed curves) are shown together with a smoothed version of the Pacific Decadal Oscillation index (red). The various curves have been mean-corrected so that they could be plotted together. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

We have mentioned several similarities in the data between the subregions such as the similarity in the decompositions obtained from EEMD. Although not shown the spectra from the various sub-regions were very similar as well. We further address the question of similarity by calculating cross-correlation functions (CCFs) between each subregion and the mean data set obtained by calculating the mean value for each month over all seven individual records. We used prewhitened versions of the original data from each location together with the mean data set (also pre-whitened) to perform the calculations. Pre-whitening removes autocorrelation from the individual time series which tends to inflate the cross-correlations. As stated earlier, this

L.C. Breaker et al. / Journal of Marine Systems 164 (2016) 151–164

161

Fig. 11. Smoothed versions of the data using LOWESS (α = 0.42) for two of the northern-most sub-regions at 1°N, 83°W (red) and 1°N, 81°W (green), with (solid), and without (dashed) the likely influence of the Pacific Decadal Oscillation (PDO) are shown in Fig. 11a. In Fig. 11b, smoothed versions of the data from the sub-regions at 1°S, 83°W (black) and 3°S, 83°W (magenta), with, and without the likely influence of the PDO are shown, using the same notation. See text for details. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

procedure is highly recommended in the analysis of climatic data (Katz, 1988). The results are summarized in Fig. 12. The pre-whitened CCFs are similar with maximum values generally falling in the range of 0.5 to 0.7, except for one departure at 1°N, 79°W. Separate calculations between individual locations confirmed that the data from this location was not well-correlated with most of its neighbors even though the results were statistically significant. This sub-region is dominated by land and because of its northerly exposure it has greater contact with, and thus influence from, the Panama Current which transports warmer waters southward along the coast of Ecuador. This sub-region has the

highest mean temperature of all the sub-regions consistent with its location and exposure (Table 1). It is noteworthy that cross-correlation analysis was far more helpful than the other methods of analysis in bringing to our attention the differences that exist between this subregion and the other six sub-regions. The mean data for all seven sub-regions is shown Fig. 13. The solid curve in red represents the long-term trend based on an EEMD decomposition of the mean data. As before, modes 8, 9, and 10 have been grouped to obtain this result. The dashed blue curve represents the corresponding result from LOWESS with α = 0.45. The dotted curves provide a measure of uncertainty and correspond to ± one standard

Fig. 12. Maximum cross-correlations for the data between each sub-region and the mean of all seven sub-regions. The data have been pre-whitened in each case (see text for details). The different colors refer to the locations of the different sub-regions. The upper 95% confidence limit has been included and is shown by the dashed line above zero. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

162

L.C. Breaker et al. / Journal of Marine Systems 164 (2016) 151–164

Fig. 13. Mean sea surface temperature for all seven sub-regions is shown by (1) an EEMD decomposition (modes 8 + 9 + 10) - the solid red curve, and (2), from LOWESS (blue) with α = 0.45. The uncertainty is shown by the dotted red curves and represents ± one standard deviation about the mean. The uncertainties were smoothed using LOWESS with α = 0.30 in order to emphasize the increasing uncertainty in the early years of the records. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

deviation about the monthly mean values, after smoothing with LOWESS (α = 0.30). We used 0.30 for the smoothing parameter in this case (less smoothing) to emphasize the increase in uncertainty at the beginning of the data where there are fewer observations. Other measures of uncertainty were considered but most provide a single value for the confidence interval and in this case a single value was clearly not appropriate. The mean curve shows that temperatures for the region increase from about 24.0 °C at the beginning to about 25.0 °C by the end of the record. We see once again the large increase in temperature between 1900 and at least 1940 – an increase that approaches 1 °C in slightly b50 years. After 1940 and up to about 1990, temperatures decrease slightly until the 1960s and then gradually increase until about 1990. From ~1990 to 2015, SSTs in the mean remain almost constant with a slight tendency to decrease by the end of the record. The influence of the first major maximum in the PDO index occurs circa 1940 in close agreement with the observed maximum in SST at approximately the same time. The decrease in SST between the 1950s and the 1980s coincides, at least approximately, with the major minimum that occurs in the PDO index from the 1950s to ~ 1980. A major minimum during this period is one of the defining characteristics of the PDO index from the last century. Finally, the change in SST from increasing to almost constant between ~ 1980 and 2000 is in phase, or consistent with, the second maximum in the PDO that occurred circa 1990.

8. Summary, discussion and conclusions Monthly-averaged SST data were acquired from seven adjacent 2° × 2° sub-regions that extend from the coast of Ecuador out to 84°W longitude, and from 2°N to 4°S latitude. The data cover the period from January 1900 through December 2014, and thus span 115 years. Based on the work of Folland and Parker (1995) a correction has been applied to the data between 1900 and 1941. This correction is related to the wooden and canvas buckets that were used to store the seawater samples that were collected prior to the actual time when the measurements of SST were taken. After 1941 the procedures were changed to avoid this problem. The correction itself represents a positive bias that increases with time reaching a value of +0.42 °C by the end of 1941. The effect of this correction has been to increase the temperatures prior to 1942 relative to those after this time which reduces

the overall long-term positive trends which we have found in the data by up to several tenths of a degree. The results of linear trend analysis indicate that the trends have positive slopes that are statistically significant in 5 out of the 7 sub-regions. In the sub-regions centered at 1°S, 83°W and 1°S, 81°W, the trends were positive but the linear slopes were not statistically significant. These trends were also used to estimate the range over which temperatures had increased at each location (Table 1). Two patterns emerged from this analysis. First, the greatest increases in SST occur in the subregions next to the coast. Second, greater increases occur at the northern-most and southern-most sub-regions while the smallest increases occurred at 1°S, 83°W and 1°S, 81°W. A second trend analysis was conducted using Ensemble Empirical Mode Decomposition (EEMD) to decompose the data into a sequence of independent modes. The three highest modes were then combined to obtain the long-term trends. By including the three highest modes, the variability at interdecadal time scales was retained. LOWESS, a state-of-the-art smoothing function, was also used to evaluate the results obtained from EEMD. The two methods produced very similar results. The results from EEMD and supported by LOWESS indicate that SST increased significantly from 1900 to at least 1940 in the majority of cases. The rate of increase in temperature during this period was increased slightly due to the application of the aforementioned correction. Also, it is likely that the Pacific Decadal Oscillation (PDO) contributed to the positive slopes between 1900 and 1940 where they occurred since the PDO and the trend were in phase during this period. Although SSTs generally increased during this time, the manner in which these increases took place differed in certain sub-regions. During the period between approximately 1950 and 1990, SSTs tend to fluctuate with slightly decreasing values followed by increasing values and a minimum that occurs during the 1960s or 1970s. After ~ 1990, both EEMD and LOWESS indicate that no major changes in surface temperature occur. The EEMD results also show that the trend for the sub-region centered at 3°S, 81°W appears to be negative unlike the other trends which all appear to be positive although weakly positive in some cases. As mentioned earlier, this sub-region is strongly influenced by the waters from the Gulf of Guayaquil which constitute a significant portion of the sub-region. The Gulf in turn receives a major influx of warmer, less saline water from the Guayas River which flows directly into the Gulf. This input from the Guayas River may be as much as

L.C. Breaker et al. / Journal of Marine Systems 164 (2016) 151–164

20 km3/year (Twilley et al., 2001). According to Stevenson (1981), the mean annual cycle of SST gradually increases from the entrance of the Gulf to the inner Gulf toward Guayaquil with temperatures generally falling in the range of 24 °C to 26 °C. Monthly maps of SST produced by the Ecuadorian Navy (INOCAR) for the coast of Ecuador typically show temperatures in the Gulf of Guayaquil to be at least 1 °C to 2 °C warmer than temperatures farther off the coast. As a result the mean temperatures for this sub-region reflect not only the influence of the oceanic waters further offshore but also of the warmer waters that typically reside inside the Gulf of Guayaquil. Thus differences in the longterm trend between this sub-region and the other sub-regions can almost certainly be attributed to major riverine influence primarily from the Guayas River. The EEMD results additionally showed that for the three sub-regions closest to the coast there is a tendency for temperatures to decrease during the first decade before they started to increase. These three subregions each have a coastal boundary which could provide insight into this behavior. It is possible that during this early period there were fewer reports from ships at sea and thus the influence of observations collected at the shoreline would influence the monthly-averaged values if they were generally higher (or lower) than those collected aboard ship. This also raises the question concerning the impact of the correction we have applied to the data? If shoreline observations are included in these data and their numbers are significant during this period, then we have the problem of often applying a correction where none is needed. Corrections during the first decade would be on the order of + 0.25 °C to + 0.30 °C. Since the correction is positive this situation would indeed increase the temperatures. Finally, it is also possible that riverine influence contributed to the warmer temperatures during the first decade at 3°S, 81°W and even in the sub-region at 1°S, 81°W, allowing for the advection of these waters further north along the coast. Clearly, this explanation, if it has merit, does not explain the anomaly that occurred at 1°N, 79°W during this period. Spectral analysis of the data was performed using two well-known methods, the periodogram and the Maximum Taper Method (MTM) of calculating power spectral densities. For these calculations we used the mean time series obtained by averaging the records from each of the 7 sub-regions. An analysis of the spectra obtained using both methods shows strong variability between 3 and 10 years, consistent with ENSO events which typically occur every 3 to 7 years. At even lower frequencies a spectral peak with a period of approximately 14 years occurs which we associate with the more recently-observed ENSO modulation cycle whose expected period is 10–15 years. A major maximum in the spectrum was also observed at approximately 0.02 cpy with a period of about 50 years or slightly greater which may well be due to the Pacific Decadal Oscillation (PDO). Although we cannot be certain that this peak is related to the PDO we are not aware of any other oscillation whose primary period falls in this range. Finally, spectral analysis of the data also revealed side band frequencies associated with the annual cycle. Side bands indicate the presence of long-period oscillations that modulate the amplitude of the annual cycle. The period associated with this oscillation is about 73 years, an oscillation that may contribute to the long-term trends we have observed. Oscillations in the range of 70 to 100 years are frequently attributed to the Gleissberg sun spot cycle and such a connection has often been referred to in the meteorological literature (Hoyt and Schatten, 1997). Focusing next on the PDO per se, although the PDO index is not formally defined at latitudes b 20°N, previous work strongly suggests that its influence is present in equatorial regions. In fact, this region may turn out to be the source of the PDO, according to Shakun and Shaman (2009). From our own work, the likely influence of the PDO was further examined using the results from Ensemble Empirical Mode Decomposition (EEMD). The 8th mode of the decomposition from each sub-region suggested the possible influence of the PDO when compared with the PDO index (Fig. 10). Then the 8th mode was subtracted from the

163

original data at each location to isolate the likely influence of the PDO on the long-term trends (Fig. 11). By subtracting mode 8 from the original data we have employed EEMD as a band-rejection filter to reduce the possible influence of the PDO. Differences between the original data and the data with mode 8 removed were generally consistent with the influence of an oscillation whose characteristics are similar to the PDO. This influence was also generally more apparent at locations further offshore. Returning to the influence of the ENSO phenomenon, it is difficult to find any obvious relationship between the trends we have observed and specific ENSO episodes. However, recent research indicates that the PDO may in large measure owe its existence to ENSO. The modeling work of Newman et al. (2003), Schneider and Cornuelle (2005), and Shakun and Shaman (2009) indicate that the PDO is primarily forced by the ENSO phenomenon. Although the models differ in detail, the bottom line is clear – ENSO is the primary driving mechanism for the PDO. Thus, ENSO, through its strong connection to the PDO, does in fact contribute in a major way to the trends we have observed. The overall similarity in the EEMD decompositions (the same number of modes in each case) and a cross-correlation analysis between the data from each sub-region suggest that the data in this region tend to be spatially homogeneous, with two possible exceptions: the sub-regions centered at 3°S, 81°W and at 1°N, 79°W. These discrepancies notwithstanding, we took the regional mean data set and likewise subjected it to an EEMD decomposition (Fig. 13). The trend extracted from the EEMD decomposition shows that temperatures for the region increase from about 24.0 °C in 1900 to approximately 25.0 °C by 2015. Once again the relatively large increase in temperature between 1900 and at least 1940 is apparent. During this period, SSTs increase at a rate approaching 1 °C in slightly b 50 years. After 1950, temperatures decrease slightly until the 1960s and then gradually increase until about 1990. From about 1990 to 2015, average temperatures remain rather constant with no significant changes during this period. The likely influence of the PDO since about 1990 may be responsible for providing sufficient cooling to offset any warming tendencies that may have occurred during this period. The likely influence of the first major maximum in the PDO index occurs circa 1940 in close agreement with the observed maximum in SST at approximately the same time. The decrease in SST between the 1950s and ~1980 coincides, at least approximately, with the major minimum that occurs in the PDO index from the 1950s to ~1980. Finally, the change in SST from increasing to almost constant between ~1980 and 2000 is in phase, or consistent with, the second maximum in the PDO that occurred circa 1990. Several conclusions can be drawn from this work and a number of questions asked. Records that are 115 years in length have provided us with the opportunity to learn about the long-term behavior in SST off the coast of Ecuador. The data indicate that the long-term variability tends to be coherent over the seven sub-regions. How dependent is this interpretation on the resolution of the data? If we had equivalent data from sub-regions that were 1° × 1° instead of 2° × 2°, would this have changed our conclusions significantly? All we can state is that our results are simply relative to the sampling schemes that were employed in creating these datasets. Although the long-term trends are not dominant features in the data, had the period of observation been limited from 1900 to 1940 or 1950, the trends would have almost certainly been dominant. The record length employed in this study has allowed us to put the last 25 years in perspective – during no period has surface temperature been as constant as it has been since ~1990. If other factors such as global warming have been at work to increase surface temperatures during this period, the cooling influence of the PDO may have served as a moderating influence. Our results depend heavily on our definition of the long-term trend. In using EEMD to obtain the trend, had we decided, for example, to include only modes 9 and 10, instead of modes 8, 9, and 10, the likely

164

L.C. Breaker et al. / Journal of Marine Systems 164 (2016) 151–164

impact of the Pacific Decadal Oscillation might have gone completely undetected. Our choice, by including slightly shorter time scales, i.e., interdecadal, has provided more insight into the oceanic processes involved, but at the expense of emphasizing more gradual movement in the data. The apparent influence of the PDO on the long-term trends cannot be overstated. The various analyses that were conducted consistently suggested the influence of the PDO, particularly during the periods where maxima or minima in the PDO index occur. Finally, although individual ENSO episodes did not clearly influence the long-term trends, recent model results indicate that they do contribute strongly over longer time scales through their contribution to the PDO. Acknowledgments The authors are especially grateful to Steve Worley and Zaihua Ji from NCAR for providing the SST data from ICOADS that has been used in this study. We also thank two anonymous reviewers for many helpful comments. Finally, we thank Dr. Allan Clarke from Florida State University for his advice and wise counsel on matters related to ENSO and tropical oceanography. References Behringer, D.W., Ji, M., Leetma, A., 1998. An improved coupled model for ENSO prediction and implications for ocean initialization. Part 1: the ocean data assimilation system. Mon. Weather Rev. 126, 1013–1021. Breaker, L.C., Loor, H., Carroll, D., 2015. Trends in sea surface temperature off the coast of Ecuador since 1900. La Tecnica 15, 100–113. Chatfield, C., 1995. The Analysis of Time Series. Texts in Statistical Science. Chapman & Hall/CRC, Boca Raton. Clarke, A.J., 1983. The reflection of equatorial waves from oceanic boundaries. J. Phys. Oceanogr. 13 (7), 1193–1207. Cleveland, W.S., 1979. Robust locally weighted regression and smoothing scatterplots. J. Am. Stat. Assoc. 74, 829–836. Cleveland, W.S., 1985. The Elements of Graphing Data. Wadsworth Advanced Books and Software, Monterey, California (323 pages). Cucalon, E., 1987. Oceanographic variability off Ecuador associated with an El Niño event in 1982-83. J. Geophys. Res. 92, 14,309–14,322. Cucalon, E., 1989. Oceanographic characteristics off the coast of Ecuador. In: Olsen, S., Arriaga, L. (Eds.), A Sustainable Shrimp Mariculture Industry for Ecuador. Narragansett, RI: Coastal Center, University of Rhode Island. Technical Report Series TR-E-6, International Coastal Resources Management Project. Davis, J.C., 1986. Statistics and Data Analysis in Geology. John Wiley & Sons, New York. Folland, K.C., Parker, D.E., 1995. Correction of instrumental biases in historical sea surface temperature data. Quart. J. Roy. Meteorol. Soc. 121, 319–367. Fuller, W.A., 1996. Introduction to Statistical Time Series, Wiley Series in Probability and Statistics. John Wiley & Sons, New York. Hoyt, D.V., Schatten, K., 1997. The Role of the Sun in Climate Change. Oxford University Press, Inc., New York.

Huang, N.E., 2005a. Introduction to Hilbert-Huang transform and some recent developments. In: Huang, N.E., Attoh-Okine, N.A. (Eds.), The Hilbert-Huang Transform in Engineering. Taylor & Francis, Boca Raton, pp. 1–24. Huang, N.E., 2005b. Introduction to the Hilbert-Huang transform and its related mathematical problems. In: Huang, N.E., Shen, S.S.P. (Eds.), The Hilbert-Huang Transform and its Applications. World Scientific, New Jersey, pp. 1–26. Huang, N.E., et al., 1998. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. London A]–>Proc. R. Soc. Lond. A 454, 903–995. Jenkins, G.M., Watts, D.G., 1968. Spectral Analysis and its Applications. Holden-Day, San Francisco. Katz, R.W., 1988. Use of cross correlations in the search for teleconnections. J. Climatol. 8, 241–253. Kessler, W.S., 2006. The Circulation of the Eastern Tropical Pacific: A Review 69 pp. 181–217. Laevastu, T., LaFond, E.C., 1970. Oceanic Fronts and their Seasonal Positions on the Surface. Naval Undersea Research and Development Center. Ocean Sciences Department, San Diego, CA. Lukas, R., 1986. The termination of the equatorial undercurrent in the eastern Pacific. Prog. Oceanogr. 16, 63–90. Mantua, N.J., Hare, S.R., Zhang, Y., Wallace, J.M., Francis, R.C., 1997. A Pacific interdecadal oscillation with impacts on salmon production. Bull. Am. Meteorol. Soc. 78, 1069–1079 (//dx.doi). Martinez-Ortiz, J., Aires de Salva, A.M., Lennert Cody, C.E., Maunder, M.N., 2015. The Ecuadorian artisanal fishery for large pelagics: species composition and spatiotemporal dynamics. PLoS ONE http://dx.doi.org/10.1371/journal.pone.0135136. Natrella, M.G., 1963. Experimental statistics. National Bureau of Standards Handbook 91. Department of Commerce, U.S. Newman, M., Gilbert, G.P., Alexander, M.A., 2003. ENSPO-forced variability of the Pacific decadal oscillation. J. Clim. 16, 3853–3857. Pak, H., Zaneveld, J.R., 1974. Equatorial front in the eastern Pacific Ocean. J. Phys. Oceanogr. 4, 570–578. Philander, S.G., 1990. El Niño, La Nina, and the Southern Oscillation. Academic Press, San Diego, CA. Sarachik, E.S., 1978. Tropical sea surface temperature: An interactive one-dimensional atmosphere-ocean model. Dyn. Atmos. Oceans 2, 455–469. Schneider, N., Cornuelle, B.D., 2005. The forcing of the Pacific decadal oscillation. J. Clim. 18, 4355–4373. Shakun, J.D., Shaman, J., 2009. Tropical origins of North and South Pacific decadal variability. Geophys. Res. Lett. 36, L19711. http://dx.doi.org/10.1029/2009GL040313. Stevenson, M.R., 1981. Seasonal variations in the Gulf of Quayaquil, a tropical estuary. Instituto Nacional de Pesca, Boletin Cientifico Y Tecnico vol. IV (No. 1). Sun, F., Yu, J.Y., 2009. A 10–15 year modulation cycle of ENSO intensity. J. Clim. 22, 1718–1735. Thompson, D.J., 1982. Spectrum estimation and harmonic analysis. IEEE Proc. 70, 1055–1096. Twilley, R.R., et al., 2001. The Gulf of Quayaquil and the Guayas River estuary. Coastal Marine Ecosystems in Latin America. In: Seeliger, U., Kjerfve, B. (Eds.), Ecological Studies vol. 144. Springer-Verlag, Berlin, pp. 245–263. Weedon, G.P., 2003. Time-Series Analysis and Cyclostratigraphy. Cambridge University Press, Cambridge, UK. Wu, Z., Huang, N.E., 2009. Ensemble empirical mode decomposition: a noise-assisted data analysis method. Adv. Adapt. Data Anal. 1, 1–41. Wyrtki, K., 1965. Surface currents of the eastern tropical Pacific Ocean. Bulletin Vol. IX (5), Inter-American Tropical Tuna Commission, La Jolla, California.