Temporal Synchrony in Satellite-Derived Ocean Parameters in the Inner Sea of Chiloé, Northern Patagonia, Chile

Abstract

Spatial synchrony occurs when geographically separated time series exhibit correlated temporal variability. Studies of synchrony between different environmental variables within marine ecosystems worldwide have highlighted the extent of system responses to exogenous large-scale forcing. However, these spatial connections remain largely unstudied in marine systems, particularly complex coastlines, where a paucity of field observations precludes the analysis of time series. Here, we used time-frequency analyses based on wavelet and wavelet coherence (WC) analysis to quantify the synchrony (co-variations) between environmental time series derived from MODIS (moderate resolution imaging spectroradiometer) in the topographically complex inner sea of Chiloé (ISC, 41–44°S) for the 2003–2022 period. We find that the strength of the synchrony between chlorophyll a ($Chla$) and turbid river plumes (for which we use remote sensing reflectance at 645 nm, $Rr{s}_{645}$) varies between the northern and southern areas of the ISC; higher synchrony, measured as the WC between these variables, is observed along the northern basin where water and particle exchanges with the Pacific Ocean are reduced. The WC analysis showed higher synchrony between these variables, with dominant periodicities of 0.5 and 1 year resulting from the hydrological regime of the freshwater input in the area that persisted throughout the 2004–2018 period. Our results suggest that the strong and significant spatial synchrony at the regional scale is likely related to the phases of large-scale climatic oscillations, as inferred through the partial wavelet coherence analysis. Potential mechanisms driving spatial synchrony are discussed in the context of climate and oceanographic regimes in the area.

1. Introduction

The study of synchrony is fundamental for understanding and describing spatial and temporal interactions in ecological and biological systems [1,2]. Spatial synchrony occurs when local populations at different geographic locations exhibit correlated dynamics over time, which can be driven by a wide variety of biotic and abiotic factors [3,4]. In addition to its spatial structure, synchrony also highlights the fundamental spatial and temporal nature of ecosystem response. Thus, time series from multiple environmental parameters can provide important insights onto the mechanisms driving synchrony [5,6]. Importantly, studies focused on the synchrony between patterns of abundance, growth, and survival rates in disjoint populations of the same species have shown that climatic oscillations are a key driver of correlated dynamics (i.e., the Moran effect) [7,8]. Previous studies of synchrony in marine ecosystems have highlighted the extent of responses to exogenous large-scale forcing [7,9]. For example, Hunt Jr et al. [10] found that climate effects on the timing of the retreat of sea ice impacted the recruitment of zooplankton on the eastern Bering Sea. Similar synchronous patterns between climatic variability and different marine groups (e.g., zooplankton, fish, and corals) have been observed over large spatial scales [11,12]. Similarly, Defriez et al. [13] argue that changes in plankton synchrony are climate-driven over a regional scale across the North Pacific. More recently, studies have identified multiple synchrony scales in environmental fields by examining patterns of spatial autocorrelation, empirical orthogonal functions (EOFs), and multiple regression techniques (e.g., [7,11,14,15,16]). In this context, the use of wavelet coherence has emerged as an efficient tool to resolve the underlying hypothesis that two signals oscillate simultaneously as is particularly well adapted to detect patterns in non-stationary time series as is common in the marine environment [15,17,18,19,20]. For example, Ménard et al. [21] explored the association between climate variability and tuna catching rates in the equatorial Indian Ocean, suggesting that climate variability and fisheries are strongly associated in periodic modes of 4–5 years. The synchrony between zooplankton abundance and climate signals at seasonal scales confirms the importance of upwelling/downwelling indices and river outflow processes in the dynamics of marine zooplankton communities [22].

The northern Patagonia region of the Southeastern Pacific ocean is a complex and highly variable environment, both in space and time (Figure 1A). Over local scales, nearshore surface water characteristics are impacted by the transition of pristine landscapes (e.g., native rainforest) to agricultural and tree plantations, changing the integrity and quality of freshwater outflows [23] and altering their relative contributions of organic matter into the coastal ocean [24]. Significant inputs of terrestrial organic carbon to the nearshore can affect the supply of mineralized nutrients, which in addition to light availability, modulate the structure of autotrophic biomass [25,26]. On the other hand, spatial gradients of surface properties impact the biogeochemical processing of inorganic nutrients, hence the temporal dynamics in primary productivity [27]. Large-scale climatic control is also important in the regulation of physical–biological coupling in the region [28,29]. For example, patterns of primary and secondary production appear to be strongly modulated by the magnitude and phase of the El Niño–southern oscillation (ENSO) cycle [30]. Moreover, owing to the exposure to low-level westerly flow, the Southern annular mode (SAM), which characterizes the north–south movement of the westerlies, has important impacts over the region [31]. During positive (negative) phases of SAM (ENSO), the equatorward intensification of the westerlies increases the advection of Subantartic waters (SAAW), a climatic pattern that has been linked to changes in the coastal oceanography, and concomitantly in ecological patterns along the inner sea of Chiloé (ISC) [28,32,33]. The spatial structure of coastal ocean circulation along northern Patagonia departs from classic meridionally oriented upwelling fronts by exhibiting marked zonally oriented fronts together with smaller-scale frontal structures [34]. The seasonal synchrony of sea surface temperature (SST) gradients with periods of high environmental variability results from seasonal river discharges and coastal upwelling events described recently by [34]. These authors identified strong SST gradients in the coastal ocean in summer and austral fall in association with a meridional coastal upwelling front. On the other hand, SST gradients are ubiquitous during spring and summer [34].

The protected waters of the ISC are of enormous importance for the Chilean aquaculture industry, which is among the top five producers globally, and accounts for almost the totality of bivalve production and a large fraction of finfish production, chiefly Atlantic salmon [35,36]. Considering the social, economic, and environmental importance of the region, establishing the factors driving synchronous dynamics in local oceanographic characteristics across the region remains a priority. However, the challenging environmental characteristics highlighted above pose a formidable challenge for the collection of physical, biological, and biogeochemical oceanographic information from the field. To date, most field studies are confined to a few protected sectors [24,25,26,37], where studies relying heavily on remotely sensed information have contributed in filling the gap [28,34]. For example, using the time series of satellite products for the marine (surface chlorophyll—a concentration, Chla; and normalized fluorescence line height, nFLH) and terrestrial (enhanced vegetation index) environments, Lara et al. [28] showed that the synchrony between these series and climatic indices (e.g., pacific decadal oscillation, PDO; El Niño–southern oscillation, ENSO) depends on the time scale, suggesting the existence of complex dynamics between phenology and climatic oscillations. The spatiotemporal fluctuations in Chla and nFLH along the meridional gradient of the ISC show higher biological activity north of the Desertores islands, a region influenced by several freshwater discharges [38,39]. The temporal synchrony between environmental time series has been less investigated due to the nonlinear nature of environmental variables and the possibility of resonance of seasonal mechanisms (e.g., SST fronts) [34] in the interaction with the complex topography. Although climate-induced changes have facilitated the inference of particular timescales of synchrony, such as the seasonal coupling between larvae supply and chlorophyll a variability [33], the synchrony dynamics between relevant ocean color parameters (as a proxy of co-variations between river discharge and surface autotrophic biomass) at a regional context have not been investigated.

This study aims to characterize the seasonal and interannual co-evolution of surface environmental characteristics along the inner sea of Chiloé (ISC) in northern Patagonia. Using satellite-derived ocean characteristics and time-frequency analysis, we take the opportunity to the hypothesis that the strength of the synchrony between environmental characteristics varies between basins of the ISC with contrasting local environmental conditions.The paper is organized as follows: Section 2 describes the data and methods, and Section 3 presents the main results. A discussion of the principal modes of synchrony and the evolution of correlation between temporal series across multiples scales is presented in Section 4. Finally, the main conclusions are presented in Section 5.

2. Materials and Methods

2.1. Study Area

The ISC is located in the northern Chilean Patagonia between Chacao Channel and Boca del Guafo (Figure 1A). It is characterized by the presence of a number of semi-enclosed basins (e.g., inlets, fjords, and estuaries) exhibiting high heterogeneity in the local environmental conditions [33,34]. The hydrography is strongly influenced by large glacial rivers (e.g., the Puelo, Futaleufú, Yelcho, and Palena Rivers, see Figure 1) that discharge considerable amounts of freshwater which maintain turbid plumes fueling biological activity in the coastal ocean, especially in summer [29,40].The lowest (highest) freshwater input to the ISC is observed during austral summer (winter–spring) [39]. In addition, this region is characterized by the confluence of the modified Subantarctic water (MSAAW, <34.0 PSU; >2 mlL⁻¹O²) and experiences intrusions of oceanic waters through Boca del Guafo at the poleward end of the ISC. The main bathymetric constriction corresponds to the Desertores islands (see Figure 1B) which separate the ISC into northern (NB) and southern (SB) basins [33,41]. On the other hand, at the west coast of Chiloé Island, two narrow submarine canyons (Ancud and Cucao) that cross the continental platform fracture the sea floor, showing depths >1000 m (see Figure 1B). The ISC bathymetric information shows a contrast between NB and SB. The first (NB) has depths > 300 m that reach maximum values along basin, while the SB is relatively shallow compared to the NB, showing depths <200 m. This meridional separation is reflected in the seasonal variability of sea surface temperature (and associated fronts) and biological activity [34,38]. In particular, the highest mean Chla and nFLH are confined to the NB (7.7 mg m⁻³ and 0.36 W m⁻² μm ⁻¹ sr⁻¹, respectively), while the lowest mean Chla and nFLH are registered in the SB where it connects with the open coastal ocean (3.8 mg m⁻³ and 0.24 W m⁻² μm ⁻¹ sr⁻¹, respectively) (Figure 1A).

2.2. Satellite Data

High spatial resolution (1000 m) daily chlorophyll a concentration (Chla; mg m⁻³), fluorescence line height (nFLH; W m⁻² μm ⁻¹ sr ⁻¹), and surface reflectance at 645 nm ($Rr{s}_{645}$) data from the resolution imaging spectroradiometer (MODIS, onboard Aqua) were used to characterize the temporal variation in environmental surface characteristics along the ISC. We used MODIS $Rr{s}_{645}$ data as a proxy for turbid (sediment-dominated) river plumes because of the low penetration in the water column and high correlation with the river discharges (e.g., [39]). The MODIS level-two data were obtained from NASA’s ocean color website (http://oceancolor.gsfc.nasa.gov, accessed on 24 May 2021). Standard procedures (e.g., atmospheric corrections) were used for Chla, nFLH, and $Rr{s}_{645}$ estimates [42]. The quality-controlled daily swaths were used to create gridded monthly composites, which were used to perform the empirical orthogonal function (EOF), wavelet, and Hilbert–Huang transform analyses (see Section 2.4).

2.3. River Discharge and Climatic Indices

Daily river discharge data were obtained from the Dirección General de Aguas de Chile (DGA) database (https://dga.mop.gob.cl/Paginas/default.aspx, accessed on 22 July 2022). We considered the stations in the Puelo and Futaleufú Rivers (Figure 1B) covering the period 2002–2020. These rivers are the greatest fluvial inputs into the northern (Puelo River) and southern (Futaleufu River) ISC basins, respectively. We considered two large-scale climate indices, the multivariate ENSO Index (MEIv2 https://psl.noaa.gov/enso/mei/, accessed on 22 July 2022) and the southern annular mode (SAM https://legacy.bas.ac.uk/met/gjma/sam.html, accessed on 22 July 2022). These climatic indices have been linked to climate and hydrographic variability along northern Patagonia [33,39,43].

2.4. Wavelet and Temporal Synchrony in Northern Patagonia

We used empirical orthogonal Function (EOF) analysis to obtain the main modes of spatiotemporal variability of the Chla and Rrs${}_{645}$ fields. We chose the singular value decomposition method to compute the EOFs to avoid the large covariance matrices resulting from the size of the satellite products [44]. We only present the first two modes as these captured most of the variance for Chla and Rrs${}_{645}$.

Wavelet analysis is used to characterize the transient dynamics of satellite-derived signals due to the non-stationary properties of ecological time series [17]. The wavelet analysis extracts time-dependent amplitude cycles and their scale relative to frequency [17,45]. The relative importance of each frequency at each time was represented in the time-frequency domain, forming the local wavelet power spectrum (WPS). Thus, wavelet analysis provides information about the periodic components of a time series and how they evolve [45]. Wavelet analysis can be also used to quantify mutual dependencies between time series and their evolution. In this study we used the Morlet wavelet, which is a complex continuous wavelet given by [45]:

$$\psi \left(t\right)={\pi }^{-1/4}exp(-i{\omega }_{0}t)exp(-t^2/2)$$

(1)

The Morlet wavelet function is the product of a sinusoidal complex $exp(-i{\omega }_{0}t)$ by a Gaussian envelope $exp(-t^2/2)$ where ${\omega }_0$ is the central angular frequency of the wavelet. The term ${\pi }^{-1/4}$ is a normalization factor to ensure unit variance. Wavelet Coherence analysis was computed to quantify statistical dependencies between two time series (i.e., the northern and southern ISC signals). wavelet coherence (WC) analysis was computed to quantify statistical dependencies between two time series (i.e., the northern and southern ISC signals). The WC is defined as the wavelet cross-spectrum normalized by the spectrum of each signal:

$$W{C}_{x,y}(f,\tau )=\frac{\Vert <W_{x,y}(f,\tau )>\Vert }{\Vert <W_x{(f,\tau )>\Vert }^{1/2}{\Vert <W_y(f,\tau )>\Vert }^{1/2}}$$

(2)

where $<>$ denotes a smoothing operator in both time and space, $W_{x,y}(f,\tau )$ is the wavelet co–spectrum of the two time series $x\left(t\right)$ and $y\left(t\right)$, and $W_x(f,\tau )$ and $W_y(f,\tau )$ are the WPS of the $x\left(t\right)$ and $y\left(t\right)$, respectively. $W{C}_{y,x}(f,\tau )$ equals one when there is a perfect linear relationship at a specific time t and frequency f between the two signals, and equals zero if both series, $x\left(t\right)$ and $y\left(t\right)$, are independent (see details in [45]). Finally, we used partial wavelet coherence (PWC), which computes the coherence between two signals after controlling for the effect of other signals. In this study the PWC analysis was performed for chlorophyll a ($Chla$) and turbidity ($Rr{s}_{645}$) in the northern and southern basins of the ISC to identify coherence (synchrony) bands between the two variables after controlling for the effects of other signals (climatic indices). We adapted the inverse PWC method for Fourier analysis [46] to wavelet analysis. This inverse PWC method is based on the spectral matrix S($f,t$), where all elements correspond to the cross-wavelet spectrum for the i and j signals, and from the $nth$ time series analyzed, $y_n($t):

$$\mathbf{S}(f,t)=\left(\begin{array}{cccc}{W}_{11}(f,t)& W_{12}(f,t)& \cdots & W_{1n}(f,t)\\ W_{21}(f,t)& W_{21}(f,t)& \cdots & W_{2n}(f,t)\\ ⋮& ⋮& \ddots & ⋮\\ W_{n1}(f,t)& W_{n2}(f,t)& \cdots & W_{nn}(f,t)\end{array}\right)$$

(3)

The PWC between $y_i\left(t\right)$ and $y_k\left(t\right)$, controlled by all other time series, is defined as:

$$PW{C}_{jk\left|\right(\setminus jk)}(f,t)={\left(\frac{\mid S^{jk}\left(f,t\right){\mid }^2}{\mid S^{jj}\left(f,t\right)\mid .\mid S^{kk}\left(f,t\right)\mid }\right)}^{1/2}$$

(4)

where $S^{jk}(f,t)$ is the $(j,k)$ element of the inverse spectral matrix ${\mathbf{S}}^{-1}\left(f\right)$ and ($\setminus jk$) denote all elements except the jth and the kth. $PW{C}_{jk\left|\right(\setminus jk)}(f,t)$ captures the frequency-specific and time localized relationship between times series j and k, by excluding the effect of all other time series. All wavelet coherence calculations were performed using MatLab software following the methods described by [45,47].

Considering the non-linear and non-stationary characteristics of ecological time series (e.g., [17,18]), we used the wavelet analysis together with a Hilbert–Huang transform (HHT), proposed by Huang et al. [48]. The HHT consists on a empirical mode decomposition (EMD) and an intrinsic mode functions (IMF) to obtain instantaneous frequencies as functions of time (see details in [49]). More information on the HHT is provided as Appendix A.

3. Results

We used the EOF analysis to characterize the spatiotemporal variability of the low-passed filtered series of the satellite-derived ocean parameters along the ISC. We only present the first two modes, which accounted for 99.9 % of the total variance in the time series for both $Chla$ and $Rr{s}_{645}$ (Figure 2 and Figure 3). The principal component analysis showed that the first mode of $Chla$ explained 80.4% of the variance, and is largely associated to the annual cycle (Figure 2A). Figure 2C shows that the concentration of $Chla$ increases in spring–summer (positive values) and decreases in fall–winter (negative values), with larger variations in the northern basin (Figure 2A,C). Typically, $Chla$ reaches its maximum annual concentrations during the summer season, with the highest concentrations in the years 2014 and 2020 (Figure 2C). The lower concentrations in the winter periods show a consistent minima, and did not show the variability observed for the summer regarding the maximum values; for example, the summer concentration of the years 2010 and 2014 showed a considerable difference in the amplitude of the first mode (Figure 2C). The spatial pattern given by the first EOF shows that the annual cycle is strongest along the western coast of the northern basin (Figure 2A). In the southern basin, the highest $Chla$ concentrations are restricted to an eastern coastal band, and the spatial distribution of the high-concentration areas are clearly influenced by the river discharges (Figure 2A). The second mode of $Chla$ accounted for 19.5% of the total variance, and highlighted a spatial pattern where the northern and southern vary out of phase; that is, an increase in concentrations along the northern basin is timed with a decrease in concentrations in the southern basin, and vice versa (Figure 2B). The temporal variability of the second mode is associated to the seasonal variability, exhibiting several peaks per year (Figure 2D).

The first two EOF modes of $Rr{s}_{645}$ are shown in Figure 3. The first mode accounted for 64.1% of the total non-seasonal variance (Figure 3C). The spatial distribution associated with the first mode highlights regions of high turbidity at the river mouths and their regions of freshwater influence along the eastern margin of the ISC (Figure 3A). The temporal variability of the first mode shows an extreme event in the winter of 2008, where the freshwater input into the ISC was considerable, particularly in the southern basin (Figure 3A,C). The second mode of variability represented a 35.8% of the variance and reproduced fluctuations on a seasonal scale (Figure 3D). The winter periods (negative values) are related to events of large freshwater discharge from the rivers along the eastern margin of the ISC and from the rivers along the southwestern region of Chiloé island (Figure 3B,D). The spatial pattern shown in Figure 3B indicates a marked turbidity contrast between coastal waters along the eastern margin and deeper waters, such that when turbidity increases at the coast, it tends to decrease towards deeper ocean waters.

The WPS for $Chla$, $nFLH$, and $Rr{s}_{645}$ in the northern and southern ISC basins are presented in Figure 4. In the northern region, the main period of variability of $Chla$ and $nFLH$ was associated with the annual cycle (Figure 4A,B). The annual period is significant during 2004–2017 for $Chla$ (Figure 4A) and throughout the entire study period for $nFLH$ (Figure 4B). In the northern region, the $Rr{s}_{645}$ WPS revealed a transient and significant one-year band between 2008 and 2010 (Figure 4C). For the southern ISC basin, the wavelet power spectrum of $Chla$ showed significant energy at the annual period between 2006 and 2015 (Figure 4D). Similarly to the northern region, we found that the $nFLH$ signal was dominated by the annual mode for the entire period (Figure 4E), while the WPS of $Rr{s}_{645}$ showed a strong transient variability over annual timescales during 2008–2012 and 2018–2019 (Figure 4F).

Our results revealed spatially synchronous processes (i.e., temporal co-evolution) at regional scales in the satellite-derived environmental variables examined. In the northern region of the ISC, the temporal relationship between the $Chla$ and $Rr{s}_{645}$ signals revealed a transient and significant synchrony for a one-year period for 2008–2015 (Figure 5A). For the years 2004–2006 and 2014–2018, these signals showed significant synchrony in the two-year periods (Figure 5A). For years 2004–2006 and 2014–2018, these signals showed significant synchrony in the two-year periods (Figure 5A). In the southern region of the ISC, a one-year synchronic period was evident for 2008–2011 (Figure 5B). In both regions, we observed the presence of intra-annual (half-year) synchronic components (Figure 5A,B). The synchrony between $Chla$ in the northern ISC and the southern ISC $Rr{s}_{645}$ signal was significant and transient in the half-year, one-year, and two-year periods between 2004 and 2015 (Figure 5C). The synchrony analysis between the southern Chla and the northern $Rr{s}_{645}$ signals revealed a marked and significant association in the one-year periodic band for 2008–2015 (Figure 5D), which is similar to the results shown in Figure 5A. We finally explored the synchrony of the $Chla$ and $Rr{s}_{645}$ signals between the northern and southern regions of the ISC. The highest synchrony in $Chla$ was observed at the annual scale (one-year periodic band), which persisted throughout the 2004–2018 period (Figure 6A). A pattern of temporal synchrony between the northern and southern ISC $Rr{s}_{645}$ signals was also observed in the one-year period band only between 2004–2012 (Figure 6B).

Further analyses using the Hilbert–Huang transform (HHT) of the satellite-derived environmental time series confirm the temporal patterns detected by the wavelet analysis. Their description, analysis, and comparisons to wavelet results are presented in Appendix A.

4. Discussion

Our results based on satellite-derived environmental parameters provide important insights into previously overlooked and significant modes of synchrony in the ISC (Figure 5). The time-frequency analysis of a suite of environmental processes provided information related to local and regional processes that may determine these synchrony patterns. In general, the results are consistent with the literature showing significant periodic bands are closely associated with seasonal SST variability, which, in turn, imposes an annual/semi-annual temperature regime [33]. The wavelet analysis and the first three iMFs (and their respective Hilbert–Huang spectra) showed the ability of these statistical techniques to extract the internal characteristics of time series and evaluate their patterns of temporal synchrony (e.g., [19,51,52,53]). On the other hand, we found that the dynamics of $Chla$ are similar to that of turbid river plumes, i.e., they share the spectral characteristics with similarities at several temporal scales evidenced in the coincident multi-peak amplitude of their marginal HHT spectra within the low-frequency band in both the northern and southern ISC regions (Figure A1 and Figure A2).

Interestingly, our results revealed large variance in both $Chla$ and $Rr{s}_{645}$ signals. The fractional contribution of the annual cycle dominates variability patterns and likely originates in the dynamics of primary productivity. Putative non-marine drivers are chiefly CDOM (colored dissolved organic matter) composition, bottom reflectance, and total suspended sediments. These elements are optically independent from phytoplankton in coastal and optically-complex inner waters (Case II waters) [54,55,56]. In this way, the remotely sensed temporal dynamics are driven by the performance of ocean color algorithms in case II waters, where water constituents can vary by orders of magnitude and may (or may not) co-vary, depending on the source and composition of the suspended and dissolved materials (e.g., [57,58]).

The spatially contrasting patterns of synchrony revealed by our analyses are primarily linked to the irregular topography and bathymetric constrictions that separate the ISC into two basins, hence having two different hydrodynamic and circulation patterns [33,37,39,41]. The southern sector experiences the inflow of oceanic waters which favors the generation of horizontal density gradients that evolve in response to wind stress and freshwater outflows [40,59]. In contrast, the northern region is strongly influenced by the seasonal influx of freshwater and suspended materials. The highest riverine outflows take place during the austral winter–spring and the lowest discharges occur during austral summer (Figure 7C,D). The temporal patterns of synchrony in our results are consistent with Vásquez et al. [38] and Flores et al. [39], who reported seasonal variability in river discharges (Puelo and Futaleufú rivers) and spatial asynchronous variability between $Chla$ and $nFLH$. Previous studies have also shown that the river discharges and autotrophic biomass cycles in northern Patagonia, measured as satellite $Chla$ and $nFLH$, are associated with large-scale climatic variability through the PDO, MEI, and SAM indices, which control precipitation across southern South America from intra-seasonal to inter-annual timescales (e.g., [33,37,39]).

Patterns of primary and secondary production in the ISC also appear to be strongly modulated by the magnitude and phase of the ENSO cycle [30]. Similarly, the SAM impacts this region directly because of its control on the meridional location of the westerlies [43]. During positive (negative) phases of SAM (ENSO), the intensification of the westerlies equatorward increases the advection of Subantarctic mode water (SAAW) into the ISC. The time-varying non-stationary trend in an environmental time series hints at the presence of an intrinsic timescale [60]. The synchrony observed between 2008 and 2018 (Figure 5) supports the hypothesis of a common external forcing mechanism influencing local drivers (i.e., climatic influence mediated by local SST). In this work, we found a multi-scale synchrony between (i) $Chla$ and turbid river plumes in the northern ISC, and (ii) $Chla$ in the southern ISC and the turbid river plumes along the northern ISC, especially between 2008 and 2016 (Figure 5). The annual mode of $Chla$ appears to be driven by turbid river plumes; the decrease in freshwater input and precipitation can thus be seen to drive a shift in key oceanographic characteristics, such as vertical stratification, along the ISC [40,61]. The non-stationary trend may then be related with the effects of climate anomalies reported for the southeast Pacific Ocean during the last 10 years [23,29,43,62]. In recent years, several studies have examined the consequences of large scale oscillations in northern Patagonia (e.g., [37,43]). Importantly, changing oceanographic conditions have been linked to ecological responses across different trophic levels within the coastal planktonic food webs [32,33,63]. Giesecke et al. [32] reported an anomalous intrusion of SAAW into the ISC during the summer which was potentially associated with changes in the intensity of SAM. This particular oceanographic and climatic scenario was also reported by Lara et al. [33]. The latter showed that a large drop in the abundance of mussel larvae was associated with changes in $Chla$ in association with large-scale climatic variability. Similar relationships between climatic conditions and anomalous freshwater discharges have been related to large-scale forcing (e.g., [39]). One of the most extreme changes in ENSO and SAM was observed during 2016–2018 (see Figure 7A,B), causing regional impacts in environmental and ecological processes [23,43]. For example, anomalous harmful algal blooms in northern Patagonia coincided with positive phases of both ENSO and SAM [23,64,65]. This region is undergoing major changes in physical conditions according with projected climate change models. The lower inputs of freshwater into the ISC during the next decades, especially during austral summer and fall, will generate environmental changes with direct impacts on the coastal ecosystems and the services they provide, such as the extensive aquaculture of salmon and mussels [66].

Our wavelet coherence analysis captured the spatial synchrony in ocean color characteristics as a proxy for phytoplankton biomass and turbid river plumes. Recent evidence suggests that the coupling between ecological processes can be observed at global scales. Revising the Moran effect, Hansen et al. [67] showed that correlated changes in mean global weather conditions associated to global warming can synchronize population dynamics or biomass/abundance patterns across vast spatial scales. In marine systems, synchrony has been shown at different spatial scales. Over local (small) scales, Cavanaugh et al. [68] showed that distances smaller than 1 km correspond to spatial synchrony scales for sea urchins and young kelps of Macrocystis pyrifera. Moreover, local environmental processes, predation, and recruitment influence are the main small-scale predictors in the synchrony of population dynamics. Over regional scales, climate-driven changes in population synchrony of several taxa of zooplankton have been shown in the North Sea [13]. Finally, at the global scale, phytoplankton synchrony has been shown to be under strong Moran control [7]. In most natural systems synchrony is attributable to a combination of biophysical drivers [22]. During the time period examined by our study, several phase–antiphase periods are evident between the climatic indices (Figure 7A,C). The opposed signals may account for the synchrony observed over several frequency bands. When we examined the coherence between $Chla$ and turbid river plumes after controlling for the effect of MEI using partial wavelet coherence, we observed significant and well-defined synchrony throughout the study period, but only for the seasonal mode at both regions (Figure 8). To estimate the PWC and control for MEI, it is required that the wavelet coherence between $Chla$ and turbid river plumes is corrected by both the wavelet coherence between $Chla$ and MEI and the wavelet coherence between turbid river plumes and MEI (Equation (4)). When comparing with Figure 5, the PWC highlights that the coherence between $Chla$ and turbid river plumes occurs essentially for the annual component, suggesting that large-scale climatic oscillations such as the MEI explain the coherence for the multi-annual components. The northern Patagonia sector is an extensive, complex, and highly dynamic area that harbors myriad human activities and harbors a unique biodiversity. Key biophysical interactions over coastal waters remain largely unstudied, and the use of satellite products across multiple spatial and temporal scales together with advanced analytical techniques can provide fundamental insights to the processes controlling local oceanographic processes.

5. Conclusions

Local environmental processes in northern Patagonia are modulated by climatic variability over the southern Pacific Ocean. To the best of our knowledge, this is the first study along northern Patagonia that (i) documents the spatial synchrony in ocean color characteristics, which we use as a proxy for phytoplankton biomass and the variability of turbid river plumes, (ii) contributes towards the long-term evaluation of spatial and temporal synchrony processes at a regional scale, and (iii) interprets results from both wavelet and HHT analysis on a dataset that is intrinsically non-stationary. The temporal dynamics in ocean color characteristics vary on different scales (semi-annual and annual timescales) between two well-defined sub-regions that exhibit distinct regimes of environmental variability. We observed that the annual mode of $Chla$ is mainly driven by the variability of $Rr{s}_{645}$, particularly for the 2008–2016 period. The continuity of long-term monitoring programs is necessary to explore the synchrony between micro basins at regional scales across the timescales of climate variability, and to understand their ecological implications for this important and unique region.