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

: 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 ﬁeld 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 ﬁnd that the strength of the synchrony between chlorophyll a ( Chla ) and turbid river plumes (for which we use remote sensing reﬂectance at 645 mn, Rrs 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 Paciﬁc 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 signiﬁcant 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.


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.

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 semienclosed 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 −1 O 2 ) 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 −3 and 0.36 W m −2 µm −1 sr −1 , 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 −3 and 0.24 W m −2 µm −1 sr −1 , respectively) (Figure 1A).

Satellite Data
High spatial resolution (1000 m) daily chlorophyll a concentration (Chla; mg m −3 ), fluorescence line height (nFLH; W m −2 µm −1 sr −1 ), and surface reflectance at 645 nm (Rrs 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 Rrs 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 Rrs 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).

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 timefrequency 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]: The Morlet wavelet function is the product of a sinusoidal complex exp(−iω 0 t) by a Gaussian envelope exp(−t 2 /2) where ω 0 is the central angular frequency of the wavelet.The term π −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: where <> denotes a smoothing operator in both time and space, W x,y ( f , τ) is the wavelet co-spectrum of the two time series x(t) and y(t), and W x ( f , τ) and W y ( f , τ) are the WPS of the x(t) and y(t), respectively.WC y,x ( f , τ) 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(t) and y(t), 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 (Rrs 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): The PWC between y i (t) and y k (t), controlled by all other time series, is defined as: where S jk ( f , t) is the (j, k) element of the inverse spectral matrix S −1 ( f ) and (\jk) denote all elements except the jth and the kth.PWC jk|(\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.

Results
We used the EOF analysis to characterize the spatiotemporal variability of the lowpassed 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 Rrs 645 (Figures 2 and 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 Rrs 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 Rrs 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 Rrs 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 Rrs 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 Rrs 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 Rrs 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 Rrs 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 Rrs 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 Rrs 645 signals was also observed in the one-year period band only between 2004-2012 (Figure 6B).The black dashed lines indicate the 95% and 90% significant areas obtained by adapted bootstrapping and the cone of influence (solid gray lines) indicates the regions where the wavelet computations are not influenced by edge effects (see [50]).The colors are coded from yellow (high coherence) to blue (low coherence).
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.  [50]).The colors are coded from yellow (high coherence) to blue (low coherence).

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/semiannual 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 (Figures A1 and A2).
Interestingly, our results revealed large variance in both Chla and Rrs 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.

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 syn-chrony 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 Rrs 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.maximum power at a period of 1.17 years (Figure A1B) for the southern region.When we compare the HHT for both regions (right panels in Figure A1), the southern ISC shows a stronger temporal evolution of iMF3, while iMF1 and iMF2 present higher values at different periods, suggesting a high intra-year dependence.The same analysis but for nFLH and Rrs 645 is presented in Figures A2 and A3, respectively.For nFLH, the iMF3 shows a strong annual period over the ISC (0.99 cycles/yr), with higher amplitude in the northern region.iMF2 shows a much broader spectrum with significantly lower periods than iMF3 (Figure A3).For the Rrs 645 signal, we observed that iMF2 and iMF3 have the highest power with peaks at 0.57 yr period (1.75 cycle/yr) and 1.04 yr period (0.96 cycle/yr), respectively, in northern ISC.Remarkably, the 1 yr period represented by iMF3 exhibits a strong decay in amplitude for the period of 2012-2014 and 2017.The southern ISC shows similar spectral distributions; the maximum frequency is roughly the annual period (1.64 cycles/yr for iFM2 and 0.85 cycles/yr for iMF3, respectively).The spectrum for iMF1 is broad, with energy at all periods, and does not show a well defined peak.Similar to the northern ISC, the iMF3 shows a marked decay in the amplitude from 2012-2014 in southern ISC.

Figure 2 .
Figure 2. The dominant patterns of the inner sea of Chiloé Chla variability over the period of 2003-2020.(A) The spatial first EOF mode and (B) second EOF modes.(C,D) show the temporal evolution of the amplitude of the first and second EOF modes, respectively.

Figure 3 .
Figure 3.The dominant patterns of the inner sea of Chiloé Rrs 645 variability over the period of 2003-2020.(A) The spatial first EOF mode and (B) second EOF modes.(C,D) shows the temporal evolution of the amplitude of the first and second EOF modes, respectively.

Figure 4 .
Figure 4. Wavelet power spectrum (WPS) showing the dominant periods (in years) of variability between 2003 to 2020 times series.(A) Chla in northern region.(B) nFLH in northern region.(C) turbid river plumes in northern region.(D) Chla in southern region.(E) nFLH in southern region.(F)turbid river plumes in southern region.The black line define the cone of influence above which computations are not influenced by edge effects.The color code for power values is graded from blue (low power) to red (high power).The black dot-dashed lines indicate the 95% and 90% significant areas obtained by adapted bootstrapping[50].

Figure 5 .
Figure 5. Wavelet coherence between (A) Chla and turbid river plumes in the northern region.(B) Chla and turbid river plumes in the southern region.(C) Chla of northern region and turbid river plumes of southern region.(D) Chla of southern region and turbid river plumes of northern region.The black dashed lines indicate the 95% and 90% significant areas obtained by adapted bootstrapping and the cone of influence (solid gray lines) indicates the regions where the wavelet computations are not influenced by edge effects (see[50]).The colors are coded from yellow (high coherence) to blue (low coherence).

Figure 6 .
Figure 6.Wavelet coherence between (A) northern and southern Chla.(B) northern and southern turbid river plumes.The black dashed lines indicate the 95% and 90% significant areas obtained by adapted bootstrapping and the cone of influence (solid gray lines) indicates the regions where the wavelet computations are not influenced by edge effects (see[50]).The colors are coded from yellow (high coherence) to blue (low coherence).

Figure 8 .
Figure 8. Partial wavelet coherence for (A) northern Chla and northern turbid river plumes, (B) southern Chla and southern turbid river plumes controlled by MEI.Black dashed lines indicate the 95% and 90% significant areas obtained by adapted bootstrapping [50] and the cone of influence (solid black lines) delimits the regions where the wavelet computations are not influenced by edge effects.The color code for synchrony values is graded from blue (low synchrony) to yellow (high synchrony).Periods of variability (in years) of the y-axis are shown in log 2 scale.

Figure A1 .
Figure A1.Left panel: The 1D Hilbert-Huang Transform of the Chla signal, representing the frequency histogram for the IMF1 (red), IMF2 (green), and IMF3 (blue).The thin horizontal line represents the main peak frequency.Right Panel: 2D Hilbert-Huang Transform.The color scale (white-blue-violet) shows the strength of energy and the solid blue line represents the main peak frequency for Chla in the (A) northern and (B) southern regions of the ISC.

Figure A2 .
Figure A2.Left panel: The 1D Hilbert-Huang transform of the Rrs 645 signal, representing the frequency histogram for the IMF1 (red), IMF2 (green), and IMF3 (blue).The thin horizontal line represents the main peak frequency.Right panel: 2D Hilbert-Huang transform.The color scale (white-blue-violet) shows the strength of energy and the solid blue line represents the main peak frequency for Rrs 645 in the (A) northern and (B) southern regions of the ISC.

Figure A3 .
Figure A3.Left panel: The 1D Hilbert-Huang transform of the nFLH signal, representing the frequency histogram for IMF2 (green) and IMF3 (blue).The thin horizontal line represents the main peak frequency.Right panel: 2D Hilbert-Huang transform.The color scale (white-blue-violet) shows the strength of energy and the solid blue line represents the main peak frequency for nFLH in the (A) northern and (B) southern regions of the ISC.

Figure A4 .
Figure A4.Time series of (A) Chl-a in northern region.(B) nFLH in northern region.(C) turbid river plumes in northern region.(D) Chl-a in southern region.(E) nFLH in southern region.(F) turbid river plumes in southern region.