Nonlinear time series analysis of coastal temperatures and El Niño–Southern Oscillation events in the eastern South Pacific

We carry out a nonlinear time series analysis motivated by dynamical systems theory to investigate the links between temperatures on the eastern South Pacific coast, influenced by the Humboldt Current System, and El Niño–Southern Oscillation (ENSO) events. To this aim, we use a set of 16 oceanic and atmospheric temperature time series from Chilean coastal stations distributed between 18 and 45^∘ S. The spectral analysis indicates periodicities that can be related to both internal and external forcing, involving not only ENSO, but also the Pacific Decadal Oscillation, the Southern Annual Mode, the Quasi-Biennial Oscillation and the lunar nodal cycle. The asymptotic neural network test for chaos based on the largest global Lyapunov exponent indicates that the temperature dynamics along the Chilean coast is not chaotic. We use local Lyapunov exponents to characterize the short-term stability of the series. Using a cross-entropy test, we find that two stations in northern Chile, one oceanic (Iquique) and one atmospheric (Arica), present a significant positive cross-dependence between local Lyapunov exponents and ENSO. Iquique is the station that presents the greater number of regional characteristics and correlates with ENSO differently from the rest. The unique large-scale study area, combined with time series from hitherto unused sources (Chilean naval records), reveals the nonlinear dynamics of climate variability in Chile.

1 Introduction

The Peru–Chile or Humboldt Current System (HCS) is a notable phenomenon of the eastern South Pacific coast. It is one of the most biologically productive eastern border oceanic currents in the world (Pauly and Christensen, 1995) as a consequence of wind-driven coastal upwellings occurring at different intensity and frequency in the southeastern Pacific (Daneri et al., 2000). Such oceanic upwelling processes also impact atmospheric temperatures (Sobarzo et al., 2007) through either local or remote effects (Hormazabal et al., 2001). Atmospheric temperatures are further affected in this region by the warm (El Niño) and cold (La Niña) phases of El Niño–Southern Oscillation (ENSO) events (Vargas et al., 2007; Enfield and Allen, 1980; Pizarro and Montecinos, 2004), and there is great interest in the possible influence of anthropogenic climate change (Sydeman et al., 2014).

A notable attempt to characterize the spatial and temporal patterns of atmospheric temperatures along the Pacific coast of South America was by Falvey and Garreaud (2009). They used in situ temperature records and satellite data for the period 1979–2006 and found that in northern and central Chile at latitudes 17–37^∘ S atmospheric temperatures were cooling by about 0.20 ^∘C per decade. They ascribed this cooling to a long-term La Niña phenomenon, which is consistent with the negative temperature trend observed in the Pacific Decadal Oscillation (PDO) for the same period. However, their analysis is limited to the period after satellite data became available in the mid-1970s: that is, just when two major El Niño events, in 1982–1983 and 1997–1998, took place. The effects of these large disturbances on temperature trends are not clear and have to be accounted for to provide reliable, unbiased conclusions. This is especially true for studies conducted in highly dynamic regions like the eastern South Pacific, where complex physical land–sea–atmosphere interactions suggest nonlinear relationships between coastal upwelling, interdecadal variability of El Niño-type events and global warming (Vargas et al., 2007).

In a context of global warming, it is important to understand the spatial and temporal patterns of atmospheric temperatures in complex systems like the eastern South Pacific region, since El Niño events have great impacts on human life (Glantz, 2022), and, in particular, on human health (Kovats et al., 2003) far beyond the eastern South Pacific (Fan et al., 2017). The nonlinear approach allows us to overcome many of the limitations of the linear framework. Its foundations were laid in the early 1980s, when deterministic chaos became a very active field of research (Bradley and Kantz, 2015), including in the geosciences (see e.g. Zeng et al., 2015; Nicolis and Nicolis, 1984; Grassberger, 1986; Pierrehumbert, 1990; Vassiliadis et al., 1991; Lorenz, 1991; Cuculeanu and Lupu, 2001; Petkov et al., 2015). Elsner and Tsonis (1992) and Tsonis et al. (1993) proposed that the atmosphere can be seen as a weakly coupled system and that the finite dimension of the attractor found by some of these aforementioned studies could correspond to the size of a subsystem. The dynamics of the atmosphere and the climate system is characterized by sensitivity to initial conditions (Kalnay, 2003): namely, a small error in observing the initial conditions is exponentially amplified and this implies the impossibility of mid- and long-term forecasting. This property was already recognized in the first developments of weather forecasting (Thompson, 1963) and was associated with the nonlinear nature of deterministic dynamical systems by Lorenz (1963).

In this work we study the dynamics that governs the climatic variability of the southeastern Pacific coast, its stabilities and instabilities. We discuss its consequences and its relationship with the HCS, the Pacific anticyclone, upwelling zones and the climatic phenomena that impact this region, thereby achieving a better understanding of the system's dynamics, which is impossible to achieve using linear methods. To accomplish this task, we perform a nonlinear analysis of a unique set of daily series of ambient coastal and sea surface temperatures, distributed between 18 and 45^∘ S and spaced every 3 to 4^∘ of latitude. We test in a mathematically rigorous fashion for the presence of chaotic dynamics through the largest Lyapunov exponent. In addition, we characterize the local stability of the systems using local Lyapunov exponents (LLEs), which quantify the state-dependent predictability over finite time horizons. Also, we study how the dynamics is affected by ENSO events by comparing LLEs and the Oceanic Niño Index in different latitudes of the Chilean coast. Finally we discuss the periodicities of the time series, the type of dynamics involved and the bioclimatic characteristics found by latitude in order to understand the role played by both internal and external forcing in the climatic variability of the eastern South Pacific.

Table 1 Daily and weekly data for oceanic stations.

^∗ The data are split into three series due to the presence of gaps.

Table 2 Daily and weekly data for atmospheric stations.

^∗ The data are split into two series due to the presence of gaps.

2 Data and methods

2.1 Data pre-processing

Sea surface temperature (SST) data were obtained from the Hydrographic and Oceanographic Service of the Chilean Navy (SHOA), which has kept records for the main ports of the country since the mid-20th century. Data on atmospheric surface temperature (AST) begin from the first half of the 20th century and were provided by the Meteorological Directorate of Chile (DMC). Since the DMC weather stations are associated with the presence of airports in major cities, they usually coincide with the main ports. For all stations, the daily data were converted to weekly average temperature time series. In the case of missing data, we applied Kalman filtering for imputation using the R package imputeTS (Moritz and Bartz-Beielstein, 2017). Tables 1 and 2 list the oceanic and atmospheric stations used in this study. The first column contains the location names that range from the northern to the southern regions of Chile and the second column indicates the time range; the third and fourth columns report the sample sizes.

Figure 1 Study area and examples of time series used in this work. (a) Sea surface temperature for stations from northern to southern Chile. (b) Study area. (c) Atmospheric surface temperature for the same region. (d) Study area relative to the Humboldt Current System (HCS). For comparative purposes, all the time series plots here use the same temperature scale and timescale.

2.1.1 Study area

Along the southeastern Pacific littoral between 18 and 45^∘ S there is a marked climatic gradient from an arid zone to a rainy Mediterranean zone. Three major divisions are identifiable along the coast of Chile (Sarricolea et al., 2017): the north at 18–30^∘ S, the centre at 30–37^∘ S and the south at 37–45^∘ S. Figure 1 shows the study area and the time series of both SST and AST used in this work.

Northern Chile is characterized by a narrow continental shelf, approximately 4–5 km wide, and an arid to semi-arid climate. Coastal upwelling is permanent in this region. Antofagasta stands out, as its orientation allows the presence of a pool of persistent hot water: a full-time El Niño (Piñones et al., 2007). La Serena, at 30^∘ S, presents a medium-sized continental shelf, approximately 20–30 km wide, where the surroundings show characteristics of persistent upwelling. Central Chile has a wider continental shelf that extends about 40 km from the coast. It is characterized by a Mediterranean climate with four well-marked seasons. This system is under the influence of a coastal upwelling driven by the wind with a strong seasonal pattern and is one of the most productive regions of the ocean (Montero et al., 2007), with more than 50 % of the fish catch of Chile and 4 % of the world catch. Its hydrographical structure is strongly affected by the runoff of several rivers, whose discharges vary seasonally with a maximum during winter off Concepción at 36^∘49^′ S. Southern Chile includes the northern limit of Chilean Patagonia, one the most extensive fjord regions in the world. The climate is rainy. It covers almost 240 000 km² with extremely complex geomorphology in one of the least densely populated areas of the country. Water circulation essentially follows a two-layered estuarine flow pattern. One superficial layer of brackish water is fed by rivers and glaciers, and another marine layer, mainly Subantarctic Surface Water (SAAW), enters predominantly through a subsurface layer (Sievers and Silva, 2008).

2.1.2 Trend estimation

One recurring issue with trend estimation is that, in most situations, the detrended series is both dependent and possibly heteroskedastic; moreover, missing data are very common and disregarding these aspects leads to invalid confidence bands. Here, we follow Friedrich et al. (2020), who solve the problem by proposing a novel autoregressive wild bootstrap scheme that does not need adjustments in the presence of missing data and is particularly suitable for climatological applications. We assume that the series X_t, $t=\mathrm{1},\mathrm{\dots },n$ admits the following decomposition:

$$\begin{array}{}\text{(1)}& \begin{array}{r}{X}_{\mathrm{t}}=m(t/n)+z_{\mathrm{t}},\\ z_{\mathrm{t}}={\mathit{\sigma }}_{\mathrm{t}}{u}_{\mathrm{t}},\\ u_{\mathrm{t}}& =\sum _{j=\mathrm{0}}^{\mathrm{\infty }}{\mathit{\psi }}_j{\mathit{\epsilon }}_{t-j},{\mathit{\psi }}_{\mathrm{0}}=\mathrm{1},\end{array}\end{array}$$

where {ε_t} is a zero-mean i.i.d. (independent and identically distributed) sequence with constant variance and finite fourth moments, m(⋅) is a smooth deterministic trend function, z_t=σ_tu_t is a stochastic process with σ_t accounting for (unconditional) heteroskedasticity, and u_t is a stationary linear process. We estimate m(⋅) by means of a local polynomial smoother (loess). For the estimation of the long-run trend the degree of the smoother has been set to 1 to mitigate boundary effects; otherwise the order is 2. The confidence bands for the nonparametric estimator of the trend have been derived using the following autoregressive wild bootstrap scheme:

1. obtain the estimate of the trend $\widehat{m}\left(t\right)$ using the bandwidth $\stackrel{\mathrm{̃}}{h}$ and derive the residuals ${\widehat{z}}_{\mathrm{t}}=X_{\mathrm{t}}-\widehat{m}\left(t\right)$;

2. derive the bootstrap errors $z_{\mathrm{t}}^{\ast }={\mathit{\xi }}_{\mathrm{t}}^{\ast }{\widehat{z}}_{\mathrm{t}}$, where ${\mathit{\xi }}_{\mathrm{t}}^{\ast }=\mathit{\gamma }{\mathit{\xi }}_{t-\mathrm{1}}^{\ast }+{\mathit{\nu }}_{\mathrm{t}}^{\ast }$, ${\mathit{\nu }}_{\mathrm{t}}^{\ast }\sim \text{i.i.d. }N(\mathrm{0},\mathrm{1}-{\mathit{\gamma }}^{\mathrm{2}})$ and $\mathit{\gamma }={\mathit{\theta }}^{\mathrm{1}/l}$, with θ=0.1 and $l=\mathrm{1.75}\sqrt[\mathrm{3}]{n}$;

3. build the bootstrap series $X_{\mathrm{t}}^{\ast }=\widehat{m}(t/n)+z_{\mathrm{t}}^{\ast }$ for $t=\mathrm{1},\mathrm{\dots },n$ and estimate the trend ${\widehat{m}}^{\ast }\left(t\right)$ upon it using the same bandwidth $\stackrel{\mathrm{̃}}{h}$ as in step 1.

The scheme provides valid confidence intervals under the assumptions as in Eq. (1) where the detrended series is dependent. For more details see Friedrich et al. (2020).

2.2 Chaos and the largest Lyapunov exponent

The largest global Lyapunov exponent quantifies the sensitivity to initial conditions and is one of the hallmarks of the presence of chaos. Assume that the series is generated by the following dynamical system in ℝ^d:

$$\begin{array}{}\text{(2)}& {\mathbf{X}}_{\mathrm{t}}=F\left({\mathbf{X}}_{t-\mathrm{1}}\right)+{\mathbf{E}}_{\mathrm{t}},{\mathbf{X}}_{\mathrm{t}}\in {\mathbb{R}}^d,\end{array}$$

where E_t is an error process in ℝ^d, which can be disregarded as far the following definition is concerned; we will comment on its role later on (see Chan and Tong, 2001, for more details). Let X₀ and ${\mathbf{X}}_{\mathrm{0}}^{\prime }$ be two close initial conditions and X_n and ${\mathbf{X}}_n^{\prime }$ their values after n time steps, respectively. Then

$$∥{\mathbf{X}}_n-{\mathbf{X}}_n^{\prime }∥\approx \exp \left(n{\mathit{\lambda }}_{\mathrm{1}}\right)∥{\mathbf{X}}_{\mathrm{0}}-{\mathbf{X}}_{\mathrm{0}}^{\prime }∥,$$

where $∥\cdot ∥\cdot$ is a suitable norm and $\mathit{\delta }=({\mathbf{X}}_{\mathrm{0}}-{\mathbf{X}}_{\mathrm{0}}^{\prime })$ is a small perturbation. In other words, as n grows, the growth factor is governed by λ₁ as other terms become negligible. Hence,

$$\begin{array}{}\text{(3)}& {\mathit{\lambda }}_{\mathrm{1}}=\underset{n\to \mathrm{\infty }}{lim}\underset{\mathit{\delta }\to \mathrm{0}}{lim}{\displaystyle \frac{\mathrm{1}}{n}}\ln \left({\displaystyle \frac{∥{\mathbf{X}}_n-{\mathbf{X}}_n^{\prime }∥}{∥{\mathbf{X}}_{\mathrm{0}}-{\mathbf{X}}_{\mathrm{0}}^{\prime }∥}}\right)\end{array}$$

is the global largest Lyapunov exponent. Note that

$${\mathbf{X}}_n-{\mathbf{X}}_n^{\prime }\approx \mathit{\delta }\cdot D{F}^{\left(n\right)}\left({\mathbf{X}}_{\mathrm{0}}\right)=\mathit{\delta }\prod _{t=\mathrm{0}}^{n-\mathrm{1}}DF\left({\mathbf{X}}_{\mathrm{t}}\right)=\mathit{\delta }\prod _{t=\mathrm{0}}^{n-\mathrm{1}}{\mathbf{J}}_{n-t},$$

where F⁽ⁿ⁾(X₀) is the n-fold iteration of the map F and J_t is the Jacobian of the map F evaluated at X_t. Hence

$$\begin{array}{}\text{(4)}& {\mathit{\lambda }}_{\mathrm{1}}=\underset{n\to \mathrm{\infty }}{lim}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}n}}\ln {\mathit{\nu }}_{\mathrm{1}}\left({\mathbf{T}}_n^{\prime }{\mathbf{T}}_n\right),\end{array}$$

where ${\mathbf{T}}_n={\prod }_{t=\mathrm{0}}^{n-\mathrm{1}}{\mathbf{J}}_{n-t}$ and ν₁(A) is the largest eigenvalue of the matrix A. λ₁ measures the average rate of divergence of nearby starting trajectories, and chaotic dynamics implies λ₁>0. The two equivalent definitions of the largest Lyapunov exponent given in Eqs. (3) and (4) are reflected in the two approaches to estimating it in finite samples. The first approach refers to Eq. (3) and is called direct since it finds close pairs of state vectors and measures the average divergence of trajectories over time (Rosenstein et al., 1993; Kantz, 1994; Giannerini and Rosa, 2001). The logarithm of such a divergence is plotted over time, and, if it is possible to identify a linear scaling region, its slope is the direct estimator of λ₁. Typically, the exercise is repeated for a range of embedding dimensions and time delays to assess the reliability of the estimate. One problem with this approach is that it is sensitive to noise: it cannot distinguish between divergence due to noise and exponential divergence due to the chaotic dynamics (Giannerini and Rosa, 2004). Moreover, there are no results on the asymptotic distribution of the estimator so that it is not possible to build a proper statistical test for the hypothesis of chaos. In Giannerini and Rosa (2001) a bootstrap solution for continuous time processes was proposed, but the intrinsic sensitivity to noise of the direct estimator remains so that we only use it to cross-check the result from the second approach, the Jacobian-based estimator, which relies upon a neural network model of the map F and its Jacobian J. The estimated largest global Lyapunov exponent ${\widehat{\mathit{\lambda }}}_{\mathrm{1}}$ is obtained by plugging the estimated quantities into Eq. (4). We use a feed-forward neural network with a single layer of “hidden units” (analogue neurons). The number of units determines the model complexity (i.e. a complicated $\widehat{F}$ that comes close to interpolating the data versus a simpler $\widehat{F}$ that smooths the data) (Ellner et al., 1992). The crucial step is the selection of the neural network model, which also provides a rigorous solution to the state-space reconstruction problem (Chan and Tong, 2001). The idea is to select the embedding parameters that provide the best fit over a grid, according to the generalized Bayesian information criterion (BIC):

$$\begin{array}{}\text{(5)}& \mathrm{BIC}=\log \left({\displaystyle \frac{\mathrm{RSS}}{n}}\right)+{\displaystyle \frac{\log \left(n\right)}{n}}[\mathrm{1}+k(d+\mathrm{2}\left)\right],\end{array}$$

where RSS is the residual sum of squares, n is the sample size, d is the embedding dimension, and k is the number of units (or nodes) of the hidden layer. The generalized BIC is consistent in that it selects the true model (if it exists) with a probability of 1 as the sample size diverges. See Appendix C for more details on the neural network approach. The modelling step allows us to filter out the effect of noise, obtain a consistent estimator for λ₁ and perform a proper statistical test for chaos (Shintani and Linton, 2004). The null hypothesis tests $H_{\mathrm{0}}:{\mathit{\lambda }}_{\mathrm{1}}\le \mathrm{0}$ against the alternative hypothesis of a positive exponent $H_{\mathrm{1}}:{\mathit{\lambda }}_{\mathrm{1}}>\mathrm{0}$. This asymptotic, one-sided test is based upon the asymptotically Gaussian standardized statistic

$$\begin{array}{}\text{(6)}& Z={\displaystyle \frac{{\widehat{\mathit{\lambda }}}_M}{\sqrt{\widehat{\mathbf{\Phi }}/M}}}\underset{M\to \mathrm{\infty }}{\overset{\mathcal{L}}{⟶}}N(\mathrm{0},\mathrm{1}),\end{array}$$

where M is the number of steps ahead and $\widehat{\mathbf{\Phi }}$ is a heteroskedasticity and autocorrelation consistent (HAC) estimator of the asymptotic variance of the M-step local Lyapunov exponent estimator ${\widehat{\mathit{\lambda }}}_M$ (see next section). The null hypothesis is rejected at level α if the value of the statistic exceeds the corresponding critical value of the standard Gaussian random variable.

2.3 Local Lyapunov exponents

A positive largest global Lyapunov exponent indicates sensitive dependence on initial conditions and is an important theoretical, global indicator of the presence of chaos. However, it is defined in the double limit of infinitesimal perturbations and infinite time steps ahead; hence it poses no practical limit to the predictability of a dynamical system (Ziehmann et al., 2000), which may typically have state-space regions of enhanced predictability, coupled with zones of instability. Nonlinear dynamical systems, be they chaotic or not, are characterized by this nonuniform predictability and one way of quantifying it is through local (finite-time) Lyapunov exponents (LLEs). In this respect, local Lyapunov exponents can be used as a time series tool that provide different complementary information with respect to global Lyapunov exponents (Bailey et al., 1997; Abarbanel et al., 1992; Smith, 1992; Wolff, 1992). We assume that the data {X_t} are from a time series generated by a nonlinear autoregressive model:

$$\begin{array}{}\text{(7)}& X_{t+\mathrm{1}}=f_d(X_{\mathrm{t}},X_{t-\mathrm{1}},\mathrm{\dots },X_{t-d+\mathrm{1}})+{\mathit{\epsilon }}_{\mathrm{t}},\end{array}$$

where X_t∈ℝ and {ε_t} is a sequence of independent random variables or perturbations with E(ε_t)=0 and Var(ε_t)=σ². It is important to note that the error in Eq. (7) is not measurement error, but dynamical noise, an inherent part of the dynamics of the system. The state-space representation of the system is

$$\begin{array}{}\text{(8)}& {\mathbf{X}}_{t+\mathrm{1}}=F\left({\mathbf{X}}_{\mathrm{t}}\right)+{\mathbf{E}}_{\mathrm{t}},\end{array}$$

where ${\mathbf{X}}_{\mathrm{t}}=(X_{\mathrm{t}},X_{t-\mathrm{1}},\mathrm{\dots },X_{t-d+\mathrm{1}})$ and ${\mathbf{E}}_{\mathrm{t}}=({\mathit{\epsilon }}_{\mathrm{t}},\mathrm{0},\mathrm{0},\mathrm{\dots },\mathrm{0})$ are in ℝ^d.

The local exponent (LLE) is defined by taking a small perturbation of X_t to ${\mathbf{X}}_{\mathrm{t}}^{\ast }$ and following forward the perturbed and unperturbed trajectories after M time steps:

$$\begin{array}{}\text{(9)}& {\mathit{\lambda }}_M\left({\mathbf{X}}_{\mathrm{t}}\right)={\displaystyle \frac{\mathrm{1}}{M}}\log {\displaystyle \frac{∥{\mathbf{X}}_{t+M}-{\mathbf{X}}_{t+M}^{\ast }∥}{∥{\mathbf{X}}_{\mathrm{t}}-{\mathbf{X}}_{\mathrm{t}}^{\ast }∥}}.\end{array}$$

If we choose to perturb only the first component of X_t and follow the growth of that perturbation along the actual sample path, this corresponds to a perturbation of Eq. (7) at time t, i.e. perturbing X_t and the following the growth of that single perturbation along the actual sample path. This definition of the LLE is related to noise amplification by the system and thus to predictability (Yao and Tong, 1994). This is best shown in a scalar system with small additive noise:

$$X_t=f\left(X_{t-\mathrm{1}}\right)+{\mathit{\epsilon }}_{\mathrm{t}}.$$

Now let ${\mathit{\sigma }}_M^{\mathrm{2}}\left(x\right)=\mathrm{Var}(X_{t+M}\mid X_{\mathrm{t}}=x)$. It results in

$$\begin{array}{c}{\displaystyle }{\mathit{\sigma }}_M^{\mathrm{2}}\left(x\right)={\mathit{\sigma }}^{\mathrm{2}}{\mathit{\mu }}_M\left(x\right)(\mathrm{1}+o(\mathrm{1}\left)\right)\text{ where }{\mathit{\sigma }}^{\mathrm{2}}=\mathrm{Var}\left({\mathit{\epsilon }}_{\mathrm{t}}\right),\text{ and}\\ {\displaystyle }{\mathit{\mu }}_M\left(x\right)=\mathrm{1}+\sum _{j=\mathrm{1}}^{M-\mathrm{1}}{\left\{\prod _{k=j}^{M-\mathrm{1}}Df\left(f^{\left(k\right)}\right(x\left)\right)\right\}}^{\mathrm{2}}.\end{array}$$

This can be rewritten as

$$\begin{array}{}\text{(10)}& {\mathit{\mu }}_M\left(x\right)=\mathrm{1}+\sum _{j=\mathrm{1}}^{M-\mathrm{1}}\exp \left\{\mathrm{2}(M-j){\mathit{\lambda }}_{M-j}\left(f^{\left(j\right)}\right(x\left)\right)\right\},\end{array}$$

where λ_j represents LLEs along the “skeleton” system with the noise deleted (Yao and Tong, 1994). However, up to the leading order of σ in the previous expansion, the skeleton LLEs agree with those along the sample path for any fixed finite M. For a system with dynamical noise, the values of μ_M set the limits to large-sample forecasting accuracy (Tong, 1995). Thus Eq. (10) shows a close link – for small noise – between LLEs and local unpredictability.

In the general d-dimensional scenario, as perturbations become small, λ_M(X_t) will depend on the derivatives of the map F. Then

$$\begin{array}{}\text{(11)}& {\mathit{\lambda }}_M\left({\mathbf{X}}_{\mathrm{t}}\right)={\displaystyle \frac{\mathrm{1}}{M}}\ln \Vert {\mathbf{J}}_{M+t-\mathrm{1}}{\mathbf{J}}_{M+t-\mathrm{2}},\mathrm{\dots },{\mathbf{J}}_{\mathrm{t}}{U}_{\mathrm{0}}\Vert ,\end{array}$$

where J_t is the Jacobian matrix of F evaluated at X_t, U₀ is a unit vector and $\Vert \cdot \Vert$ is some vector norm. Since λ_M(X_t) is a function of the state-space location, the LLE depends on the trajectory and can be thought of as an M-step-ahead local Lyapunov exponent process. If U₀ is chosen at random with respect to the uniform measure on the unit sphere, then with a probability of 1 as M→∞, λ_M(X_t) converges to the global Lyapunov exponent λ₁ because U₀ has zero probability of falling into a subspace corresponding to subdominant exponents. In practice, we usually take $U_{\mathrm{0}}=(\mathrm{1},\mathrm{0},\mathrm{0},\mathrm{\dots },\mathrm{0})$ and the estimators for local exponents are the same as for global ones, with the notable difference that the former depend upon the state-space location and measure the state-dependent instability of the system.

Figure 2 Spectral analysis for three representative stations from the northern, central and southern regions of Chile. Column (a) corresponds to sea surface temperature time series (oceanic stations) and column (b) to atmospheric surface temperature time series (atmospheric stations). Dominant peaks are indicated together with their associated periods in years.

3 Results

We carry out a univariate spectral analysis for the 16 stations. In Figure 2, we show the results for three atmospheric and three oceanic stations, which represent the northern, central and southern regions of Chile. The complete set of results is in Appendix B. The oceanic stations that show spectral peaks with the longest periodicities are Valparaíso with 27 years, Talcahuano with 20 years, Antofagasta with 18 years and Arica with 14 years. Four stations have periodicities with an intermediate value of 6 to 10.6 years. Three stations (Antofagasta, Valparaíso and Talcahuano) present a 4-year periodicity. Two stations have peaks around 3 years, and six stations present 2-year periods. In the eight stations, periodicities between 1.3 and 2 years were also found: most stations show periods of around 1.5 years.

The atmospheric stations that show spectral peaks with the longest periodicities are Antofagasta and Serena with 20 years, Arica with 19 years, Iquique with 18.8 years, Concepción with 16.4 years, and Puerto Montt with 16 years. Furthermore, Arica and Antofagasta share a periodicity of 5.1 years. There are three stations with periodicity around 4 years and five stations with periodicity around 3 years. Periodicities around 1.5 years were also found, as in the oceanic case. In general, longer periodicities can be observed for atmospheric surface temperature time series than for sea surface temperature. Possible causes of these periodicities can be found in the discussion in Sect. 4.1.

Figure 3 Temperature time series of Antofagasta and Iquique. (a) Antofagasta sea surface temperature; (b) Antofagasta atmospheric surface temperature. (c) Iquique sea surface temperature; (d) Iquique atmospheric surface temperature. Interpolated original (light blue), trend cycle (blue) and long-term trend (black), together with confidence bands at 95 % (light blue).

We estimate the trend by means of the local polynomial smoother as detailed in Sect. 2.1.2. For the estimation of the long-term trend the degree of the polynomial was set to 1 to mitigate boundary effects; otherwise the order is 2. The confidence bands at the 95 % level for the nonparametric estimator of the trend were derived using an autoregressive wild bootstrap scheme, which is robust against the presence of missing data, serial dependence and heteroskedasticity. We show the results for Antofagasta (SST, AST) and Iquique (SST, AST) in Fig. 3. The estimated long-term trend for Antofagasta and Iquique might seem to indicate a tendency, but this is ruled out if we take into account the uncertainty (confidence bands in light blue). This suggests that the process that has generated the time series for both stations is mean-stationary, and this holds for all the stations. For more details, see Appendix B1.

Table 3 Output from the best neural network fit from the grid search test for chaos of the time series of the oceanic stations. The table reports the embedding dimension d, the time delay τ, the number of hidden units k, the value of the minimized BIC, the global Lyapunov exponent ${\widehat{\mathit{\lambda }}}_{\mathrm{1}}$ (together with the value of the test statistic Z^∗) and the p value of the test.

Table 4 Output from the best neural network fit from the grid search test for chaos of the time series of the atmospheric stations. The table reports the embedding dimension d, the time delay τ, the number of hidden units k, the value of the minimized BIC, the global Lyapunov exponent ${\widehat{\mathit{\lambda }}}_{\mathrm{1}}$ (together with the value of the test statistic Z^∗) and the p value of the test.

We estimate the largest Lyapunov exponent of the series and derive a rigorous (asymptotic) statistical test for chaos by following Shintani and Linton (2004). The approach is based upon a single-layer, feed-forward neural network model of the map and its Jacobian. More details on the procedure are reported in Sect. 2.2 and Appendix C. The null hypothesis of no chaos ($H_{\mathrm{0}}:{\mathit{\lambda }}_{\mathrm{1}}\le \mathrm{0}$) is tested against the alternative hypothesis of a positive exponent ($H_{\mathrm{1}}:{\mathit{\lambda }}_{\mathrm{1}}>\mathrm{0}$). One of the critical issues is the selection of the best neural network model, which we solve by means of the generalized Bayesian information criterion (BIC), where the global minimum is searched on a grid of embedding dimensions from d=3 to d=20 and time delays from τ=1 to τ=4. The results are presented in Tables 3 and 4 for the oceanic and atmospheric stations, respectively. Besides the station name and the time period, the tables report the values of the embedding dimension d, time delay τ and number of hidden units k of the network that produce the minimized BIC criterion over the parameter grid. The last two columns contain the value of the test statistic Z^∗ and the p value of the test based on the asymptotic Gaussian null distribution, respectively. All the estimated exponents are negative and the test never rejects the null hypothesis of the absence of chaos for both the oceanic and the atmospheric stations.

Figure 4 Attractors for sea surface temperature with superimposed information on the local Lyapunov exponents. Red corresponds to a positive LLE, whereas green indicates a negative LLE. The subsections indicate the latitudinal order of the stations. In the case of Arica and Talcahuano there are three courts, which are numbered from 1 to 3. The temporal windows are 1 month for Iquique; 4 months for Arica; 5 months for Coquimbo, Valparaíso and Talcahuano; and 6 months for Antofagasta, Corral and Puerto Montt.

Figure 5 Attractors for atmospheric surface temperature with superimposed information on the local Lyapunov exponents. Red corresponds to a positive LLE, whereas green indicates a negative LLE. The subsections indicate the latitudinal order of the stations. For the case of Rodelillo, there are two courts. The temporal windows are 5 months for Serena, Rodelillo, Concepción and Valdivia and 6 months for Arica, Iquique, Antofagasta and Puerto Montt.

This result is an important indication that the dynamics behind temperature in Chile is most probably nonlinear but not chaotic (Lorenz, 1963). We use local Lyapunov exponents (LLEs) to characterize local predictability throughout the stations. We varied the number of steps ahead, M, from 4 to 52 (i.e. from 1 month to 1 year) for oceanic stations and up to 260 steps ahead (5 years) for atmospheric stations. In order to uncover any association with El Niño, we chose three specific windows corresponding to the short (0 to 3 months), medium (6 to 9 months) and long term (12 to 22 months) and paid particular attention to the results obtained for the 1-, 6- and 12-month windows. In Figs. D1 to D3 of the Appendix we show the time series of the oceanic stations where the information on the 1-month, 6-month and 1-year LLEs, respectively, is superimposed: red corresponds to a positive LLE, whereas green indicates a negative LLE. The same information for the atmospheric stations is reported in Figs. D4 to D6. Likewise, in Figs. 4 and 5, we present the phase-space portraits of the attractors for the stations with superimposed information on the LLEs in red and green.

Figure 6 Box plots of the local Lyapunov exponents ${\widehat{\mathit{\lambda }}}_M$ versus steps ahead (M in weeks) for the sea surface temperature series. The series are arranged latitudinally.

Figure 7 Box plots of the local Lyapunov exponents ${\widehat{\mathit{\lambda }}}_M$ versus steps ahead (M in weeks) for the atmospheric surface temperature series. The series are arranged latitudinally.

Comparing the LLEs of the oceanic stations with their atmospheric counterparts for the different windows, we found that for the 1-month window, both classes of stations showed unstable behaviour throughout the record, except for Iquique SST that shows patterns of both stability and instability. For the 6-month window in SST and AST, mixed stability and instability is observed, except for Iquique, where there is stability throughout the record. In the Antofagasta and Puerto Montt SSTs similar behaviour can be observed, showing instability at the extremes and stability in the middle part. Finally, for the 1-year window, all oceanic stations showed stability throughout the record, except for Talcahuano and Puerto Montt: that is, the southern region. Of the atmospheric stations, Arica, Iquique and Antofagasta showed instability in the extremes and stability in the middle part of the record; Rodelillo behaves totally stably, while Serena, Concepción and Valdivia present totally stable or unstable parts throughout the series. This is not the case for Puerto Montt, which presents very little instability. The box plots of the LLEs versus steps ahead (M) are shown in Figs. 6 and 7. Clearly, the sea surface temperature series approach the global exponent faster than the atmospheric series. For instance, for M=28 steps ahead, the box plots of the LLEs are already negative for all the SST series, with the exception of Puerto Montt (Fig. 6). On the contrary, for M=28 all AST time series still have positive LLEs (Fig. 7). Overall, the results indicate a qualitative difference in the temperature dynamics of the two types of series. This can be explained by considering that, on average, ${\widehat{\mathit{\lambda }}}_{\mathrm{1}}\approx -\mathrm{0.05}$ for the sea surface series but ${\widehat{\mathit{\lambda }}}_{\mathrm{1}}\approx -\mathrm{0.01}$ for the atmospheric series, so the latter are globally more stable. The presence of the ocean and its exchanges of heat and momentum with the atmosphere can easily reduce the instability properties of the atmospheric flow (Vannitsem et al., 2015).

In the following, we compare the information conveyed by local Lyapunov exponents with the Oceanic Niño Index (ONI). El Niño is characterized by a positive ONI greater than or equal to +0.5 ^∘C. La Niña is characterized by a negative ONI less than or equal to −0.5 ^∘C. To be classified as a fully fledged El Niño or La Niña episode, these thresholds must be exceeded for a period of at least five consecutive overlapping 3-month seasons. The Climate Prediction Center of the US National Oceanic and Atmospheric Administration (NOAA) considers El Niño or La Niña conditions to occur when the monthly Niño 3.4 SST departures meet or exceed ±0.5 ^∘C along with consistent atmospheric features. These anomalies must also be forecast and should persist for 3 consecutive months. In the case of the oceanic stations, we note the following: in Antofagasta the expected behaviour is observed – that is, the whole series shows persistent instability for the 1- and 6-month windows throughout the record (Fig. D1c). This might be due to the warm water pool that is always found here simulating an El Niño. The same pattern is observed for the 1-year window, but, in this case, the behaviour is consistently stable (Fig. D3c). Arica shows instability throughout the record for the 1-month window, where instability decreases during El Niño events of 1953–1954, 1958–1959, 1963–1964, 1965–1966, 2009–2010 and 2014–2016 (Fig. D1a). During the 6-month and 1-year windows, it becomes more stable than unstable and shares the same behaviour with Talcahuano, where the third cut can be located lower than the other two and a very marked variability in the LLE can be seen (Figs. D2a, f and D3a, f). At these two stations the trajectory of the phase space for the second cut is interrupted (Fig. 4). In Talcahuano, interesting behaviour, with respect to the ONI, is only observed for the 1-month window, where the instability for the first cut for 1949–1974 decreases during El Niño events of 1953–1954, 1958–1959, 1963–1964, 1965–1966, 1968–1969 and 1972–1973, to show rising behaviour in the third cut of 1991–2020 during El Niño events of 1991–1992, 1994–1995, 1997–1998; 2002–2003, 2004–2005, 2009–2010 and 2015–2016 (Fig. D1f). In the case of Iquique, the 1-month and 1-year windows have stable and unstable behaviour; it is the only series in which the peaks of both the stable and unstable part coincide with several El Niño events of 1987–1988, 1991–1992, 1994–1995, 1997–1998, 2002–2003, 2004–2005, 2006–2007, 2009–2010, 2014–2016, 2018–2019 and 2019–2020 (Figs. D1b and D3b). For the 6-month window, the behaviour is stable throughout the record, but the stability peaks continue to coincide with the same El Niño events, with the exception of 2002–2003, 2004–2005, 2006–2007, 2018–2019 and 2019–2020 (Fig. D2b). Finally, Valparaíso, in the 1-month window, shows a decrease in the unstable part that corresponds to El Niño events of 1968–1969, 1972–1973, 1982–1983, 1997–1998 and 2014–2016 (Fig. D1e).

As concerns atmospheric stations, in Antofagasta the 1-month window shows unstable behaviour throughout the entire record and shows a temperature decrease during the 1972–1973, 1982–1983 and 1997–1998 El Niño events (see the highest temperature peaks in Fig. D4c). In the case of Arica, for the same window, a decrease in the unstable part is observed for the same El Niño events plus that of 2014–2016 (see the highest temperature peaks in Fig. D4a). Puerto Montt behaves in the same way as Antofagasta in its oceanic record: that is, it shows continuous unstable behaviour throughout the entire time span (Figs. D1h and D4h). The same occurs with Rodelillo for the 1-month window (Fig. D4e); its behaviour is similar to that of Antofagasta and Puerto Montt. Finally, for the 1-year window in Serena, instability increases from 1971 to 1982–1983, which is where it coincides with the El Niño event of that year. The behaviour decreases and increases again to coincide with El Niño of 1987–1988, and likewise, a peak of instability coincides with El Niño of 2014–2016 (Fig. D6d).

In summary, in the case of the oceanic stations, in Arica, Valparaíso and Talcahuano, the unstable part decreases significantly during various El Niño events, while for Antofagasta the behaviour is persistent throughout the record for both the stable and the unstable part; on the other hand, in Iquique the unstable behaviour increases during El Niño events and the stability decreases after them. For the atmospheric stations, it is in Arica, Iquique and Antofagasta where a decrease in instability during several El Niño events can be seen, contrary to what Serena shows with an increase in instability; only the latter shows this for the 1-year window, while the previous stations show an interesting behaviour of decrease for the 1-month window.

Figure 8 Cross-entropy S_k for $k=-\mathrm{12},\mathrm{\dots },\mathrm{12}$ between the 1-month and 6-month LLEs (Iquique SST) and between the 1-month LLEs (Arica AST) and the ONI. (a, c, e) Original time series (b, d, f) pre-whitened time series. Iquique 1 month and 6 months correspond to the first two upper graphs and the two central ones, while Arica 1 month corresponds to the lower two graphs.

In order to go beyond this visual comparison and to corroborate whether there is a relationship between the ONI and the dynamics of climate variability in Chile seen through the LLEs, we compute the cross-entropy measure proposed in Giannerini et al. (2015) and Giannerini and Goracci (2023); see Appendix D. Rejection bands under the null hypothesis of cross-independence are derived through the moving block bootstrap, and the series have been pre-whitened using an autoregressive filter in order avoid spurious rejections due to serial dependence. The test was performed for all stations before and after pre-whitening. The most significant results are shown in Fig. 8, whereas the full set of results can be found in Appendix E. For Iquique SST, the cross-entropy test indicates a significant dependence at lag +1 and −2 with 1-month and 6-month LLEs, respectively. For the atmospheric stations, only Arica AST showed a significant dependence at lags +2 and −3. This implies that the ONI at time t is correlated with the LLE at times t+1, t−2, t+2 and t−3, respectively, for the aforementioned stations. For the 6-month time horizon, Corral also presents a significant dependence, which could be a false positive; see Fig. E4. We also applied the test to correlate the ONI and 1-year-ahead LLEs for Serena. The results in Fig. E7 do not show any significant cross-dependence. We offer a further discussion in Sect. 4.2.3 below.

4 Discussion

4.1 Spectral analysis

The periods of 18, 18.8 and 19 years found for Antofagasta SST, Iquique AST and Arica AST could be linked to the PDO. They might also be associated with the lunar nodal cycle of 18.61 years, popularly known as “lunar wobble” – this low-frequency tidal cycle is due to the oscillation of the maximum lunar declination with respect to the Equator (Rossiter, 1962; Black et al., 2009) – since the variations of the lunar cycle influence the movement of ocean currents. The 4-year periods found in Antofagasta, Valparaíso and Talcahuano might be linked to the lunar perigee subharmonic of 4.4 years. The external forcing caused by the Moon is important, not only because of the changes caused in interannual timescales to the variability of the tidal range, but also because it influences the hydrodynamic energy of the ocean, altering sedimentation patterns (Oost et al., 1993; Gratiot et al., 2008) and coastal erosion (Smith et al., 2007), water level and coastal flood risk (Thompson et al., 2021), and ocean stratification (Devlin et al., 2017b, a, 2019).

Interannual periods can also be observed in the oceanic stations: 2.7 years for Iquique and Valparaíso, 2.6 for Arica, 2.5 for Coquimbo and Talcahuano, 2.4 for Antofagasta, 2.0 for Corral, and 1.8 for Puerto Montt. In the case of atmospheric stations, 2.7 years was found for Serena, 2.8 for Rodelillo, 2.3 for Valdivia and 2.0 for Puerto Montt. These periods could be related to the Quasi-Biennial Oscillation (QBO), a climatic alteration located in the Pacific Ocean that consists of zonally symmetrical easterly and westerly wind regimes that regularly alternate with a period that varies between 24 and 30 months and is manifested through hydroclimatological variables such as wind and temperature. Its eastern phases appear to be alternately related to the warm component of ENSO (Riveros et al., 2020). The QBO has already been seen in other records from the Chilean region between 39 and 43^∘ S (Villalba et al., 1998). Fagel et al. (2008), when analysing layered sediments from Lake Peyuhue (40^∘ S, 72^∘ W), found that the frequency bands of the thickness of the layers were mainly related to QBO periodicities. In our case, we found this oscillation in oceanic and atmospheric surface temperature records. For the oceanic records it was found in all the stations, something that did not happen with the atmospheric records, where it is only seen for the stations in the central and southern parts of Chile. Similarly, periods related to ENSO and PDO were obtained, as well as some others that we believe may be related to local events; however, it should be noted that QBO has also been related to ENSO.

Regarding the periodicity of 1.5 years found in both classes of station, Jiang et al. (2020) have noted this periodicity in modes of tropical Pacific surface wind anomalies, also referred to as the so-called combination mode (C-mode). This C-mode emerges from nonlinear interactions between the warm pool annual cycle and ENSO variability (Stuecker et al., 2013). Sakova et al. (2006) likewise found an 18-month signal in altimeter sea surface height over the whole Indian Ocean, which they associate with Indian Ocean dipole events, similar to ENSO, and statistically linked to it (Stuecker et al., 2017), but taking place in the equatorial Indian Ocean.

4.2 Nonlinear analysis and ONI

The effects of changes in the SST and AST on biota range from the redistribution of species (Lubchenco et al., 1993; Ihlow et al., 2015), the increase in vulnerability that occurs in the HCS due to the expansion of oxygen minima associated with the increase in SST (Mayol et al., 2012; Stramma et al., 2008), and increases in extreme events (Harris et al., 2018; Garreaud, 2018) to subtle changes in biological–climatic feedbacks (Bradford et al., 2019; Bardgett et al., 2008), all of which invite discussions about the contribution of the oceans to these feedbacks (Joos, 2015). For these reasons, understanding the effects of climate variabilities such as El Niño on the behaviour of long time series of coastal temperature records can contribute to improving forecasts and understanding the effects of climate change and climate variability in systems such as the HCS (Harris et al., 2018). The results consistently indicate the presence of nonlinear, non-chaotic oscillators for both oceanic and atmospheric stations. LLE analysis highlights patterns of reduced and enhanced stability and predictability in the phase space, where the asymptotic behaviour of the LLEs towards the negative global Lyapunov exponent differs in oceanic versus atmospheric stations, since stability is reached faster in oceanic stations, where the oceanic environment is less perturbable than the atmosphere due to the great heat capacity of water.

4.2.1 Oceanic stations

For both oceanic and atmospheric stations, Antofagasta shows the slowest speed of convergence in the oceanic case, accompanied by Corral. Antofagasta is characterized by the constant presence of warm water inside its bay (Piñones et al., 2007), so its climate possesses a type of local permanent El Niño, keeping the bay in a state of external forcing. This makes the dynamics of the Antofagasta system look unstable in the LLEs, and when comparing with ONI it is not possible to distinguish different behaviour when a general El Niño event occurs. Regarding Corral, there is a cooling trend in the coastal upwelling regime along northern Chile and southern Peru, which contrasts with the warming trend in the last 350 years in an offshore area of the Humboldt Current System at 36^∘ S (Vargas et al., 2007). Such El Niño-type conditions also found in Corral might explain its reduced speed of convergence. The case of Iquique SST is peculiar not only due to its attractor, which shows an abrupt change, but also because of the estimated embedding dimension (d=3) that is much lower than that of the other stations, with mixed dynamics for the window of 1 month and a positive cross-dependence between the ONI and the LLE. This indicates a strong relationship between the climate variability of this region and the El Niño phenomenon. Fuenzalida (1985), through an analysis of El Niño 1982–1983, showed that an increase in temperature in Iquique is associated with this phenomenon. In our case, the cross-dependence test showed a significant correlation at lag +1 for the 1-month LLEs and at lag −2 for the 6-month LLEs. The difference between the two is that the record observed at 1 month presents stabilities and instabilities, which, together with other factors discussed below, favour a leading effect and could be associated with the warm phase of ENSO. As for the 6-month LLEs, the behaviour of the area is completely stable; perhaps for this reason there is a lagging effect associated with the cold phase of ENSO. Notably, like La Serena and Valdivia, with similar dynamics occurring between them as discussed below, Iquique is also in a transition zone, the faunal transition zone, located south of the subtropical and north of the subantarctic zone (Bé and Tolderlund, 1971). In addition, the northern part of Chile, in which Iquique is located, is an area with wind-driven coastal upwelling throughout the year. Barbieri et al. (1995) analysed the northern region and found four sources of upwelling in Iquique. The upwelling south of Iquique is the most important due to its latitudinal extension. For this area the highest correlation between wind and SST is observed. Due to the peculiar oceanic and environmental characteristics of the Iquique area, we can hypothesize that the abrupt change observed in the trajectory of the SST transition could be due to the climate of 1976–1977, which, according to Namias (1978), Houghton et al. (2001) and Solomon et al. (2007), was a coupled global variation of the atmosphere–ocean system that links the variability of the sea surface temperature of the Pacific and climatic parameters in most of the world. Although this was a global variation, this abrupt trajectory change is present in Iquique, in contrast to the other stations, where the shape of the attractors indicates the presence of nonlinear oscillators, possibly forced. Likely, these peculiar features associate Iquique's dynamics with ENSO events. Also, we found similarity between the dynamics of Arica  – upwelling throughout the year – and Talcahuano  – seasonal upwelling. This could be the result of large-scale oceanic processes that occur simultaneously in both places. In Castro et al. (2020), similar variations of δ¹³C of particulate organic matter were noted in Iquique and Talcahuano, despite being two totally different areas; this was hypothesized to be influenced by coastal outcrop. In our case, since Arica is not very far from Iquique, we can think of this outcrop as part of the explanation for the similarity of the dynamics that we found. On the other hand, we have the northward currents that carry subantarctic water, as well as the intrusion of subtropical waters to the south, which during El Niño conditions can approach the coastal zone, causing the propagation of trapped waves to the coast and a poleward subsurface flow. From a dynamical point of view, this variability can affect currents, the mixing of the water column, the intensity of upwelling, and the temperature and sea level on the continental shelf and slope throughout the Humboldt Current System (Shaffer et al., 1997, 1999). In addition to this, the Chilean continental shelf is extremely narrow compared to the Atlantic coast, with a maximum width close to 45 km in the Talcahuano area and a general maximum depth of 150 m, except for the Valparaíso area, where it reaches up to 800 m (Camus, 2001). So, these may be other possible explanations for the similarity in dynamics found in these two regions. However, something else must be happening between the two stations to maintain this dynamics in the three study windows.

4.2.2 Atmospheric stations

Rodelillo, Puerto Montt and Antofagasta ASTs present similarities with Antofagasta SST. The former three stations present persistent behaviour throughout the time span for both northern and southern regions of Chile. Therefore, we can infer the presence of teleconnection processes throughout the zone. Gómez (2008) made hydrographic and meteorological observations during the period 1996–2005, finding an interannual and decadal modulation of the beginning of the upwelling season during the period 1968–2005 for the central–southern zone of Chile. This modulation would be linked to changes in atmospheric circulation patterns, associated with teleconnections triggered during El Niño. According to this study, the presence of teleconnections in oceanographic processes is needed. Finally, for the atmospheric stations of Arica, Iquique and Antofagasta, a decrease in unstable behaviour can be seen during El Niño events, possibly due to the effects of the Pacific anticyclone. Regarding Arica, the cross-entropy test showed a significant cross-dependence between ENSO and 1-month LLEs at lags +2 and −3. Puerto Montt is the southernmost station on the northern border of Patagonia (Garreaud et al., 2013). This station is located well inside the continent, on the coast of an inland sea, with greater thermal amplitudes due to the continental effect. As this station is located south of the southern limit recognized for the influence of El Niño, it is possible that El Niño is not responsible for the climate variability at this station, but rather the Southern Annual Mode (SAM) (González-Reyes, 2013), since we can also see unstable dynamics in the oceanic part (Puerto Montt SST). However, for this region there is already a high degree of disturbance, so what could be happening here are further intensifications induced by El Niño. This may be the reason why a behaviour similar to that of Antofagasta SST is also observed when we compare the ONI and the LLE AST of both stations.

In the case of Rodelillo, we can hypothesize that the local climate is greatly influencing the results obtained; this series was divided into two parts, which showed homogeneous results for the three different types of windows: that is, if the first segment of Rodelillo for the 1-year window was stable, the second was, too. This area is in the central zone of Chile and most of the wildfire activity in the country (Soto et al., 2015; Urrutia-Jalabert et al., 2018) is concentrated in both the central and central–southern zones of Chile (32–43^∘ S), where ENSO and the Antarctic Oscillation (AAO) have a significant influence (Urrutia-Jalabert et al., 2018). This may be what maintains homogeneous results in Rodelillo's dynamics.

For their part, La Serena and Valdivia show similar dynamics for the 6-month windows. La Serena, in a transition zone between the hyperarid northern Atacama Desert and the more mesic Mediterranean region of southern central Chile, is strongly affected by ENSO, which is the primary driver underlying the high interannual variability in precipitation, and which promotes higher average rainfall during the warm phase (Rutllant and Fuenzalida, 1991). On the other hand, Valdivia is located within a transitional area between the influence of ENSO and SAM. The dynamics of both seasons is influenced by ENSO; however, the effect of the tropical South Pacific anticyclone also plays an important role in this similarity, as the effects of this semi-permanent high-pressure system change from 35^∘ S, 90^∘ W in January to 25^∘ S, 90^∘ W in July (Kalthoff et al., 2002; Montecinos et al., 2015): that is, in a period of 6 months, which is one of the windows that we use and in which important changes in the dynamics of the systems can be seen. The influence of the anticyclone is permanent to the north of 30^∘ S, maintaining dry conditions, and weak south of 40^∘ S (Camus, 2001), where there is rain throughout the year, which is why it has important effects in Valdivia. For La Serena, the cross-entropy test did not show significance; Fig. E7. This indicates that, although the instability of the dynamics in Serena is influenced by ENSO, it is not this phenomenon that contributes the most to the instabilities. We found similar results for Valdivia. This was expected because this area is also influenced by SAM and it is likely that the latter has a greater influence on the area.

4.2.3 ONI

Of all the oceanic stations tested, only Iquique SST (1- and 6-month LLEs) showed a significant cross-dependence with ENSO. In the case of atmospheric stations, only Arica AST (1-month LLEs) showed a significant cross-dependence with ENSO. This indicates that the association between LLEs and ENSO only occurs in the northern part of Chile. A possible explanation for this may be the disturbances coming from the Equator that only influence the northern part of Chile and that can affect the physical environment (Pizarro et al., 1994). On the other hand, the correlations obtained with a forward or backward lag could be related to the initial conditions that influence the westerly winds, since, according to Fedorov (2002) and Fedorov et al. (2003), during the beginning of a warm event, a westerly wind can accelerate the development of the event, while one after the peak of El Niño will simply prolong its duration. The majority of the published studies on ENSO in Chile report findings for the northern and central regions (Gutiérrez and Meserve, 2003; Holmgren et al., 2001; Guera and Portflitt, 1991; Vilina et al., 2002; Ulloa et al., 2001; Acosta-Jamett et al., 2016); although ENSO is a phenomenon that affects the climate globally, it has a greater influence on the northern part of our study region. Each El Niño event it is likely to affect the various areas of the Chilean region differently. This can be seen in Barbieri et al. (1995), who observed that during the 1987 El Niño the mean SST showed values above 1.5 ^∘C throughout the year, while during the 1992 El Niño these anomalies occurred from March to May and from September to November. It is a fact that El Niño influences the regional climate, as do the PDO and the SAM.

5 Conclusions

We examined with an LLE analysis the patterns of local instability for the 16 stations, and these correlate on some occasions with ENSO events. However, it is only for the northern stations, specifically Iquique SST and Arica AST, that a significant correlation can be observed. ENSO, despite being a phenomenon that affects the climate globally, has greater influence in the northern part of Chile.

Teleconnections are noted throughout the entire region that may be associated with ENSO, PDO, SAM, the Pacific anticyclone, upwelling zones, QBO and the C-mode. A link can be seen between the climate variability of the region and different internal forcings (PDO, ENSO, SAM, C-mode and Pacific anticyclone) and external ones (QBO, lunar cycle) that contribute to the complexity of the system and favour a change in the variability as we move in latitude, which is due, in part, to the fact that it is necessary to take regional characteristics into account if one wants to understand the response of the different studied systems, oceanic and atmospheric, to the different forcings.

Previously unused Chilean naval temperature data from a large latitudinal spread of stations have helped us to uncover a great variety of factors involved in the temperature dynamics of the Chilean region. We find a non-chaotic climate variability, but with a nonlinear behaviour, reinforcing what has been mentioned in the literature. Our study invites more detailed work on the northern part of Chile, especially on Iquique and its sea surface temperature, where the behaviour of the system was very peculiar in both its low embedding dimension and the shape of the attractor as well as the significant correlation with ENSO for time horizons of 1 and 6 months.

Appendix A: The Humboldt Current System, El Niño–Southern Oscillation and Southern Annual Mode

The Humboldt Current System (HCS) that brings water from the tropical and sub-polar regions is characterized by evident interannual and seasonal variabilities (Mayol et al., 2012). The system is composed of (i) the Humboldt Current predominantly away from the equatorial coast that moves at an average surface velocity of 6 cm s⁻¹, (ii) coastal currents formed by the Peru–Chile current and countercurrent, and (iii) the Peru–Chile equatorial coastal current. The currents towards the pole are responsible for transporting equatorial and subtropical subsurface waters to the Chilean coasts, while the flow towards the Equator brings cold Antarctic and Antarctic intermediate waters. The HCS is controlled to a large extent by the equatorial coastal winds linked to the Pacific subtropical anticyclone, which promotes the coastal upwelling associated with the primary productivity in the north and centre of Chile, extending its influence to the south of Chile in summer (Montecino et al., 2006). Around the western tropical zone of the Pacific there are ocean–atmosphere interactions that have effects at the planetary level (Enfield and Allen, 1980). In particular there is the El Niño–Southern Oscillation (ENSO) phenomenon, characterized by two opposite phases: one of warming and rainfall in the eastern Pacific, known as El Niño, and a second of cooling and dry years called La Niña (Pizarro and Montecinos, 2004). The redistribution of temperature in the water column can be associated with the passage of low-frequency coastal trapped waves (CTWs) (days at the intra-seasonal scale) and cool and warm interannual anomalies produced by equatorial waves due to the Southern Oscillation of El Niño (Shaffer et al., 1997; Montecinos and Gomez, 2010). In addition to ENSO, there is the Pacific Decadal Oscillation (PDO), which develops in the northern Pacific Ocean. It may be described as a long-period ENSO (Núñez et al., 2013), since it also has a positive or warm phase and a negative or cold phase but a different timescale. The PDO can intensify or diminish the impact of ENSO, depending on the relative phase in which these two oscillations are found. When ENOS and PDO are in phase, the impacts induced by El Niño or La Niña will be magnified with respect to normal patterns. Conversely, if ENSO and PDO are in opposite phases, the effects on global climate variability will be weaker (Wang et al., 2014). The Southern Annual Mode (SAM) or Antarctic Oscillation defines the changes in the westerly winds that are driven by atmospheric pressure contrasts, which in turn generate pressure differences between the tropics and the southern polar areas. The change in position of the western wind band, produced from west to east in latitudes between 30 and 60^∘ of both hemispheres, influences the strength and position of cold fronts and midlatitude storm systems. In the positive phases of SAM, a belt of strong westerly winds contracts towards Antarctica. This results in weaker than normal westerly winds and high pressures in southern Australia, restricting the entry of colds fronts. SAM differentially affects the surface temperature of the four continental masses of the Southern Hemisphere, including part of the area influenced by the HCS at midlatitudes (Gillett et al., 2006). In conjunction with ENSO, SAM also partially explains the coastal temperature anomalies at high latitudes (Fogt and Bromwich, 2006) and also at midlatitudes (Yeo and Kim, 1992).

Appendix B: Spectral analysis

Figures B1 and B2 show spectral analyses performed for the oceanic and atmospheric stations. Numbers are given for the spectral peaks with the highest energies.

Figure B1 Spectral analyses for oceanic stations. The numbers indicate the spectral peaks with the highest energy.

Figure B2 Spectral analyses for atmospheric stations. The numbers indicate the spectral peaks with the highest energy.

B1 Trend estimation

In Figures B3 and B4 trend estimations are given for the 16 stations in a latitudinal manner from north to south.

Figure B3 Trend estimation for oceanic stations. Interpolated original (light blue), trend cycle (blue) and long-term trend (black), together with confidence bands at 95 % (light blue).

Figure B4 Trend estimation for atmospheric stations. Interpolated original (light blue), trend cycle (blue) and long-term trend (black), together with confidence bands at 95 % (light blue).

Appendix C: Neural network models for random dynamical systems

As in Sect. 2.2, the dynamical system has the following representation:

$$\begin{array}{}\text{(C1)}& {\mathbf{X}}_{t+\mathrm{1}}=F\left({\mathbf{X}}_t\right)+{\mathbf{E}}_{t+\mathrm{1}},{\mathbf{X}}_{\mathrm{t}}\in {\mathbb{R}}^d,\end{array}$$

where E_t is an error process in ℝ^d. We assume that the data {X_t} are generated by the nonlinear autoregressive model

$$\begin{array}{}\text{(C2)}& X_{t+\mathrm{1}}=f_d(X_{\mathrm{t}},X_{t-\mathrm{1}},\mathrm{\dots },X_{t-d+\mathrm{1}})+{\mathit{\epsilon }}_{\mathrm{t}},\end{array}$$

where X_t∈ℝ and {ε_t} is a sequence of independent random variables with E(ε_t)=0 and Var(ε_t)=σ². Also, we denote with J_t the Jacobian of the map, evaluated at X_t. The model of Eq. (C1) can be seen as the state-space representation of the system of Eq. (C2), where ${\mathbf{X}}_{\mathrm{t}}=(X_{\mathrm{t}},X_{t-\mathrm{1}},\mathrm{\dots },X_{t-d+\mathrm{1}})$ and ${\mathbf{E}}_{\mathrm{t}}=({\mathit{\epsilon }}_{\mathrm{t}},\mathrm{0},\mathrm{0},\mathrm{\dots },\mathrm{0})$. We derive a consistent estimator for the map F and its Jacobian through a neural network estimator of f_d. Neural networks are a class of nonlinear models inspired by the neural architecture of the brain (Nychka et al., 1992). These are made up of layers that in turn have connected “neurons”, which send messages and share information with each other. Layers are classified into three groups: (1) input, (2) hidden and (3) output. The input values X_t are received by the input units, which simply pass the input forward to the hidden units u_j. Each connection (indicated by an arrow) performs a linear transformation determined by the connection strength ω_ij so that the total input to unit u_j is

$$\begin{array}{}\text{(C3)}& \sum _{i=\mathrm{1}}^d{\mathit{\omega }}_{ij}{X}_{t-i+\mathrm{1}}+{\mathit{\omega }}_{\mathrm{0}j},\end{array}$$

and each unit performs a nonlinear transformation on its total input:

$$\begin{array}{}\text{(C4)}& u_j=\mathit{\psi }\left(\sum _{i=\mathrm{1}}^d{\mathit{\omega }}_{ij}{X}_{t-i+\mathrm{1}}+{\mathit{\omega }}_{\mathrm{0}j}\right).\end{array}$$

The activation function ψ is a sigmoidal function with limiting values 0 and 1 as x→ −∞ and + ∞, respectively. Here we use

$$\begin{array}{}\text{(C5)}& \mathit{\psi }\left(x\right)={\displaystyle \frac{x(\mathrm{1}+|\frac{x}{\mathrm{2}}\left|\right)}{\mathrm{2}+\left|x\right|+\frac{x^{\mathrm{2}}}{\mathrm{2}}}}.\end{array}$$

Other choices of the activation function are possible even if it does not have a major impact upon the final output. The hidden layer outputs u_j are passed along to the single output unit, which performs an affine transformation on its total input. Therefore, the network output f_d, for d inputs and k units in the hidden layer, can be represented as

$$\begin{array}{}\text{(C6)}& \begin{array}{rl}{f}_d\left({\mathbf{X}}_{\mathrm{t}}\right)=& f_d(X_{\mathrm{t}},X_{t-\mathrm{1}},\mathrm{\dots },X_{t-d+\mathrm{1}})\\ =& {\mathit{\beta }}_{\mathrm{0}}+\sum _{j=\mathrm{1}}^k{\mathit{\beta }}_j{u}_j\\ =& {\mathit{\beta }}_{\mathrm{0}}+\sum _{j=\mathrm{1}}^k{\mathit{\beta }}_j\mathit{\psi }\left(\sum _{i=\mathrm{1}}^d{\mathit{\omega }}_{ij}{X}_{t-i+\mathrm{1}}+{\mathit{\omega }}_{\mathrm{0}j}\right).\end{array}\end{array}$$

We made use of the R nnet package to implement the estimator, and we used least-squares minimization (Venables and Ripley, 2002). We also experimented with conditional maximum likelihood without noticing major discrepancies in the results. The major challenge in regression with neural network models is selecting among the many possible combinations of d and k. Here we choose the best model that minimizes the BIC defined as

$$\begin{array}{}\text{(C7)}& \mathrm{BIC}=\log \left({\widehat{\mathit{\sigma }}}^{\mathrm{2}}\right)+{\displaystyle \frac{\log \left(n\right)}{n}}[\mathrm{1}+k(d+\mathrm{2}\left)\right],\end{array}$$

where the error variance is estimated through the residual sum of squares ${\widehat{\mathit{\sigma }}}^{\mathrm{2}}=n^{-\mathrm{1}}\mathrm{RSS}=n^{-\mathrm{1}}{\sum }_{t=\mathrm{1}}^n(X_t-{\widehat{f}}_d({\mathbf{X}}_{t-\mathrm{1}}){)}^{\mathrm{2}}$. Once the estimator ${\widehat{f}}_d$ is obtained, a consistent estimator for the map F and its Jacobian J can be derived with plug-in methods, as described in Shintani and Linton (2004).

Appendix D: Local Lyapunov exponents

In Figures D1 to D6 local Lyapunov exponents are shown for all 16 stations for 1-month, 6-month and 1-year windows.

Figure D1 Temperature record for ocean stations with superimposed information on 1-month local Lyapunov exponents: positive values for the LLEs are shown in red, and negative values are shown in green. Stations are given in a latitudinal arrangement.

Figure D2 Temperature record for ocean stations with superimposed information on 6-month local Lyapunov exponents: positive values for the LLEs are shown in red, and negative values are shown in green. Stations are given in a latitudinal arrangement.

Figure D3 Temperature record for ocean stations with superimposed information on 1-year local Lyapunov exponents: positive values for the LLEs are shown in red, and negative values are shown in green. Stations are given in a latitudinal arrangement.

Figure D4 Temperature record for atmospheric stations with superimposed information on 1-month local Lyapunov exponents: positive values for the LLEs are shown in red, and negative values are shown in green. Stations are given in a latitudinal arrangement.

Figure D5 Temperature record for atmospheric stations with superimposed information on 6-month local Lyapunov exponents: positive values for the LLEs are shown in red, and negative values are shown in green. Stations are given in a latitudinal arrangement.

Figure D6 Temperature record for atmospheric stations with superimposed information on 1-year local Lyapunov exponents: positive values for the LLEs are shown in red, and negative values are shown in green. Stations are given in a latitudinal arrangement.

Appendix E: Cross-entropy tests with the measure S_k

Let {X_t} and {Y_t} (t∈ℕ) be two stationary random processes, where $F_{X_{\mathrm{t}},Y_{\mathrm{t}}}(x,y)=P(X_{\mathrm{t}}\le x,Y_{\mathrm{t}}\le y)$, $F_{X_{\mathrm{t}}}\left(x\right)=P(X_{\mathrm{t}}\le x)$ and $F_{Y_{\mathrm{t}}}\left(y\right)=P(Y_{\mathrm{t}}\le y)$. Then, the metric entropy S_ρ at lag k is a normalized version of the Bhattacharya–Hellinger–Matusita distance, defined as

$$\begin{array}{}\text{(E1)}& \begin{array}{rl}{S}_k=& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\int \int \left(\sqrt{d{F}_{(X_t,Y_{t+k})}(x,y)}\right.\\ & {\left.-\sqrt{d{F}_{X_t}\left(x\right)d{F}_{Y_{t+k}}\left(y\right)}\right)}^{\mathrm{2}}\\ =& \mathrm{1}-\int \int \sqrt{d{F}_{(X_t,Y_{t+k})}(x,y)d{F}_{X_t}\left(x\right)d{F}_{Y_{t+k}}\left(y\right)}.\end{array}\end{array}$$

In the case where Y_t=X_t for all t, S_k measures the serial dependence of {X_t} at lag k and can be interpreted as a nonlinear autocorrelation and cross-correlation function that overcomes the limits of Pearson's correlation coefficient. As pointed out in Maasoumi (1993), Granger et al. (2004) and Giannerini et al. (2015), S_ρ satisfies many desirable properties, including the seven Rényi axioms and the additional properties described in Maasoumi (1993). We compute the cross-entropy tests for dependence as detailed in Giannerini and Goracci (2023). The estimator of S_k is based on univariate and bivariate nonparametric kernel density estimation and numerical integration. The confidence bands are derived under the null hypothesis of mutual independence between series. In Figures E1–E6 we show the results of the cross-entropy test between the ONI and LLEs for the oceanic stations with a time horizon of 1 and 6 months, as well as for the atmospheric stations with a time horizon of 1 month. For Serena we also tested the 1-year-ahead LLEs; see Fig. E7. For the atmospheric stations, only 1-month LLEs are shown, since no relationship was found between the ONI and 6-month-ahead LLEs.

Figure E1 Cross-entropy S_k for $k=-\mathrm{12},\mathrm{\dots },\mathrm{12}$ between the 1-month LLEs (SST) and the ONI for the first four stations. Left: original time series. Right: pre-whitened time series.

Figure E2 Cross-entropy S_k for $k=-\mathrm{12},\mathrm{\dots },\mathrm{12}$ between the 1-month LLEs (SST) and the ONI or the last four stations. Left: original time series. Right: pre-whitened time series.

Figure E3 Cross-entropy S_k for $k=-\mathrm{12},\mathrm{\dots },\mathrm{12}$ between the 6-month LLEs (SST) and the ONI for the first four stations. Left: original time series. Right: pre-whitened time series.

Figure E4 Cross-entropy S_k for $k=-\mathrm{12},\mathrm{\dots },\mathrm{12}$ between the 6-month LLEs (SST) and the ONI for the last four stations. Left: original time series. Right: pre-whitened time series.

Figure E5 Cross-entropy S_k for $k=-\mathrm{12},\mathrm{\dots },\mathrm{12}$ between the 1-month LLEs (AST) and the ONI for the first four stations. Left: original time series. Right: pre-whitened time series.

Figure E6 Cross-entropy S_k for $k=-\mathrm{12},\mathrm{\dots },\mathrm{12}$ between the 1-month LLEs (AST) and the ONI for the last four stations. Left: original time series. Right: pre-whitened time series.

Figure E7 Cross-entropy S_k for $k=-\mathrm{12},\mathrm{\dots },\mathrm{12}$ between the 1-year LLEs for La Serena and the ONI. (a) Original time series. (b) Pre-whitened time series.