Insight into anthropogenic forcing on coastal upwelling off south-central Chile

Eastern ocean boundaries are some of the most productive areas on Earth as a result of wind-driven coastal upwelling systems (CUS), which supply nutrient-rich subsurface waters to the surface layer. Upwelling nutrients support high levels of chlorophyll and high primary production rates (Chavez and Smith, 1995) which in turn support economic activities such as fisheries (Pauly and Christensen, 1995). Large-scale tropospheric subsidence that maintains the Southeast Pacific Anticyclone (SPA) influences wind intensity and direction, thus regulating Chilean CUS (Figure 1a). The effects of climate change are likely to manifest in coastal areas in various ways, and changes in CUS represent one important example of these impacts (e.g., Echevin et al., 2012; Sydeman et al., 2014; Wang et al., 2015).

Diverse areas of upwelling off the Chilean coast have been well identified and described (Strub et al., 1998; Thiel et al., 2007). Subsurface waters, particularly of equatorial origin and known as Equatorial Subsurface Water, lift cold, nutrient-rich and oxygen-poor waters to the surface. Equatorward winds over the coastal ocean are responsible for a majority of CUS by triggering Ekman pumping and offshore Ekman transport (Pickett and Paduan, 2003; Bravo et al., 2016).

CUS are also affected by wind intensity, bottom topography, coastal geometry, and other oceanographic conditions, partially associated with remote forcing, such as coastal-trapped waves or El Niño Southern Oscillation (ENSO) events (Shaffer et al., 1997, 1999). South-central Chile (35–42°S) exhibits exceptionally strong seasonal fluctuations in wind stress, precipitation and river runoff. During the austral fall–winter, the SPA moves northward, allowing for the arrival of extratropical cyclones and increasing the frequency and intensity of precipitation, freshwater discharge, and poleward winds related to coastal downwelling episodes (e.g., Sobarzo et al., 2007). Conversely, the SPA migrates southward during spring and summer, reducing precipitation and river discharge, but intensifying northward upwelling-favorable winds (e.g., Rahn and Garreaud, 2014). During this same period, more frequent and intense atmospheric low-level coastal jets can be observed on a synoptic scale, owing to an increase in the alongshore sea level pressure gradient produced by mid-latitude perturbations (Garreaud and Muñoz, 2005). As a result, the seasonal and synoptic variability of coastal winds is associated with changes in the position and intensity of the SPA, as well as sea level pressure (SLP) anomalies to the south (Muñoz and Garreaud, 2005).

Changing patterns in wind stress, precipitation and river discharge could modulate water column dynamics in terms of mixing and/or stratification (Sobarzo et al., 2007). These processes control the extent of nutrient fertilization, which modulates biological production and oxygen consumption (Escribano and Morales, 2012). In this way, coastal hypoxia appears to intensify as a consequence of climate change (Stramma et al., 2008), and it is possible that an acute dissolved oxygen deficit causes fish mortality, modifies microbial communities, C and N cycling, and stimulates the production of greenhouse gases, such as methane and nitrous oxide (Grantham et al., 2004; Naqvi et al., 2010; Gilly et al., 2013).

Since the late 1970s, an ongoing decline in precipitation has been registered in central Chile (CR2 2015; Boisier et al., submitted). This pattern has been partially attributed to changes in large-scale atmospheric circulation and Pacific Decadal Oscillation (PDO). However, it is unlikely that the entirety of this variation arises from natural causes; instead, a part of it is most likely due to the regional effects of anthropogenic climate change (Boisier et al., 2016). Furthermore, the reduction in precipitation over the Southeast Pacific is a robust signature of climate model projections for future scenarios, as well as for recent decades (Seth and Rojas, 2003). Lower precipitation rates have resulted in water scarcity and reductions in river runoff (Garreaud et al., 2017), which may result in decreased nutrient supply to coastal waters and therefore affect coastal primary productivity (Massotti I, personal communication).

There is ongoing debate regarding climate change impacts on coastal upwelling. Model simulations suggest that changes in wind stress by the end of the 21st Century as a result of climate change will increase upwelling at high latitudes in the CUS (Rykaczewski et al., 2015; Wang et al., 2015). However, whether a climate change signal can already be observed in wind stress patterns is unclear. Some regional observations support Bakun’s hypothesis (Bakun, 1990) related to land warming and ocean cooling triggered by the intensification of alongshore winds on the CUS (Santos et al., 2005, 2012; García-Reyes and Largier, 2010; Sydeman et al., 2014). Therefore, it is of primary importance to maintain reliable records and gain a better understanding of how the coastal ocean responds to anthropogenic forcing. The novel contribution of this study is to combine different approaches to better elucidate how signals of anthropogenic climate change have forced trends in upwelling-favorable winds through the spatial configuration of SLP trends over recent decades. The results allow us to explore changes in upwelling-favorable winds and their influence on the physical variables of sea surface temperature (SST) and salinity, which have been shown previously in CUS (Gutierrez et al., 2011; Aravena et al., 2014; Schneider et al., 2017), and to examine changes over the last decades in the biogeochemical variables of chlorophyll a (chl a), nutrients and dissolved oxygen concentrations off the shore of south-central Chile (30–42°S).

The atmospheric data analyzed includes monthly SLP and wind at a height of 10 m above sea level, provided by ERA-Interim reanalysis (Dee et al., 2011) and Global Climate Model (GCM) simulations from 23 models (Table S1) associated with the Coupled Model Intercomparison Project Phase 5 (CMIP5). Fully atmosphere-ocean coupled simulations and SST-forced simulations were used for 1979 to 2005. The fully coupled ensemble employs historical runs (hereafter referred to as HIST), while the SST-forced simulations emerge from the Atmospheric Model Intercomparison Project (AMIP); the latter were carried out using the atmospheric component of each of the 23 GCMs included in HIST. HIST and AMIP simulations therefore employed exactly the same climate forcing, with the exception of SST conditions. We decided to compare the reanalyzed trends with the simulated trends over the same time period, rather than using the entire historical simulations, in order to take advantage of the observational portion of the reanalyzed data; in particular, the ERA-Interim reanalysis has updated physics, improved assimilation data, and higher resolution (Dee et al., 2011), as well as better homogeneity through time, making it more reliable for modeling long-term processes (Stopa and Cheung, 2011). In order to determine if changes in CUS are attributable to climate change in fully coupled simulations (HIST), internal variability was assumed to be negligible by averaging the entire dataset to assess trends of anthropogenic climate change forcing over previous decades (Boisier et al., 2016). Simple empirical attribution was achieved by calculating the SLP trend linearly, congruently with the PDO for the AMIP simulations. First, the regression coefficient was calculated for standardized monthly SLP time series and the monthly PDO index (Mantua et al., 1997). The resulting regression coefficients were then multiplied by climate index trends, resulting in SLP trends that were linearly congruent with natural decadal forcing.

Offshore Ekman transport was used as a coastal upwelling index:

where τ_y is the alongshore wind stress, f is the Coriolis parameter and ρ_w is seawater density (1,027 kg m^–3). To derive upwelling-favorable wind stress, winds were first projected onto the alongshore direction and then the standard bulk aerodynamic formula from wind speed at a height of 10 m, using a constant drag coefficient of 0.0013 (Kraus, 1972; Nelson, 1977). We used this constant value of drag coefficient following Large and Pond (1981), who demonstrated that it does not generate large variations in the calculation of stresses when the wind speeds have values between 4 m s^–1 and 11 m s^–1. In our study area, the mean wind magnitude fluctuates within that range (Sobarzo et al., 2007). Accumulated Ekman transport from September to February was calculated to examine interannual upwelling variability.

Monthly chl a composites and photosynthetically active radiation (PAR) data, collected by the Aqua MODIS satellite for 2003–2016 with a spatial resolution of 4 km, were analyzed. In order to estimate phytoplankton biomass as chl a, 27 1° × 0.5° pixel quadrants were created between 36°S and 40°S. The average value of PAR and chl a was calculated for each time interval in the time series within each quadrant; an average estimate was considered valid if it contained at least 33% accurate pixels.

In-situ SST data recorded using near surface thermistors at four different locations in south-central Chile (Figure 1b) was provided by the Hydrographic and Oceanographic Service of the Chilean Navy (SHOA) for Coquimbo (30°S–71.4°W), Valparaíso (33°S–71.6°W), Talcahuano (36.7°S–73.1°W) and Valdivia (39.9°S–73.4°W). In addition, a fixed time-series station, known as Station 18 (St18), has been sampled at approximately 30-day intervals since October 2002, making it a well documented site (e.g., Escribano and Morales, 2012). The station is located on the continental shelf, 10 nautical miles (18.5 km) off the coast of central Chile (36.5°S; 73.1°W) at ~92 m depth (Farías et al., 2015). This time-series station is thus located in the middle of one of the most extensive continental shelves off central Chile, where it is subjected to a coastal upwelling process during austral spring and summer (Sobarzo et al., 2007). Upwelling causes strong horizontal density gradients at the sea surface (i.e., an upwelling front), limiting the active upwelling area and acting as a boundary between coastal and adjacent oceanic surface water (Letelier et al., 2009). As St18 is within the active upwelling area, it captures coastal processes occurring on the continental shelf.

In this study, continuous temperature, salinity, oxygen, fluorescence, and PAR profiles were obtained using data from Sea-Bird (SBE-19 and SBE-25) CTDs, which were reviewed and pre-processed using software from Sea-Bird Electronics, Inc. Data were processed at the original sampling frequency and averaged for every 1 dbar, then post-processed using different quality control protocols in order to minimize potential biases or errors.

Seawater at St18 was sampled using a 10-L Niskin bottle mounted on a 10-bottle rosette. Samples were obtained over ten depths, equally distributed between the surface and bottom layers, to test the concentration of gases (i.e., oxygen and nitrous oxide), nutrients, and pigments (in this correlative temporal order). Triplicate samples for nutrients (NO₃^– and PO₄^3–) were collected by connecting syringes directly to the Niskin bottle spigots and filtering through a 0.45-μm UptiDisc® adapted to the syringe. Samples were then stored at –20°C and analyzed using standard manual (from 2002 to 2009) or automated (from 2009 to 2016, SEAL Analytical AutoAnalyzer) colorimetric methods (Grasshoff et al., 1983). For chl a, 250 mL of seawater were filtered (in triplicate) using a glass-fiber filter (0.7 GF/F) which was immediately frozen until later fluorometric analysis (Turner Design AU-10), according to standard procedures (Parsons et al., 1984). For more details, see Farías et al. (2015).

Nutrient flux towards the surface (F) was calculated as follows:

where UI is the mean Upwelling Index (m³ s^–1), calculated over 100 m of coastline (Bakun, 1973), averaged over the last three days of each field campaign in the Concepción area; and C is the observed mean nutrient concentration (mmol m^–3) at depths between 40 and 80 m (subsurface layer) at St18. This simple model assumes total replacement of the zonal transport by vertical advection (i.e., without considering meridional advection).

A threshold method was implemented to determine mixed layer depth (MLD). In this method, MLD is located where temperature and density profiles change by 0.5°C and 0.15 kg m^–3 (Monterey and Levitus, 1997), relative to surface reference values for 10 and 0 m depths, respectively. This method was selected because it allows for addressing different processes that may be occurring at the same time, particularly in a dynamic coastal zone that experiences seasonal upwelling/downwelling and direct input of freshwater through rain and riverine discharge, as at St18. The criteria used provide an accurate representation of the MLD within the study zone, where the temperature criterion (T_crt) shows the variability associated with wind stress and heat fluxes, and the density criterion (D_crt) incorporates freshwater fluxes.

Monthly temperature, salinity, oxygen, nitrous oxide, chl a, PAR and nutrient concentration (both NO₃^– and PO₄^3–) anomalies were calculated as follows:

where xN is the discrete value at a certain time and XN is the mean monthly climatological value obtained using the times series between 2003 and 2016.

Time-series anomaly regressions for each variable were calculated to identify trends in oceanographic variables at St18 for the period 2002–2016, using the generalized least squares method. Trends were observed for summer (December to February) and winter (May to July), as well as for individual calendar months. To test the statistical significance of these trends, we used the non-parametric Mann-Kendall test for significance with a p-value less than 0.05. We also performed linear regression and non-parametric Spearman tests to identify correlations between variables and investigate oceanographic and biogeochemical relationships.

Natural interdecadal and interannual variability, related to the PDO and ENSO, have been demonstrated to be key contributors to the observed atmospheric and oceanographic changes off the coast of south-central Chile (e.g., Shaffer et al., 1999; Rahn, 2012). Specifically, the intensity and position of the SPA exhibit an interdecadal oscillation capable of modifying large-scale atmospheric circulation, with consequences over the CUS (Ancapichún and Garcés-Vargas, 2015). Notably, over past decades (1979–2016) the PDO shifted from a positive to a negative phase, resulting in a negative trend of the PDO index (Figure 2a). It is therefore expected that trends in SLP and winds over recent decades may be influenced by the negative PDO tendency. Nevertheless, whether changes observed off the south-central coast of Chile accurately reflect anthropogenic influences on the climate remains unclear.

Alongshore wind variability off the Chilean coast is primarily controlled by nearby SLP differences; the cross-shore SLP gradient forces geostrophic alongshore winds, while an alongshore SLP gradient produces a semi-geostrophic alongshore flow due to the presence of coastal mountains. In order to assess wind-driven changes in upwelling within the region, we first focused on the evolution of large-scale atmospheric circulation and calculated the linear trend of SLP between 1979 and 2005 (Figure 2). A clear increase of SLP is evident in the Southeast Pacific (Figure 2b). Furthermore, a spatially coherent pattern of SLP trends is observed for all data sets, although GCM simulations (HIST and AMIP) exhibit lower trends than ERA-Interim data. Major changes are observed by the SPA and near the south-central Chilean coast (36–45°S). However, the response to changes in PDO is more significant toward the south-central Pacific (Figure 2e), where it represents a slight contribution of the PDO to trends observed in SLP for the coastal region. In this case, a notably similar spatial configuration is observed using both AMIP and HIST simulations.

Assuming that SLP responds to climate change and PDO linearly, the sum of these two trends (HIST and AMIP-PDO) should be closer to the total trends. Hence, under linearity, the stronger trends observed toward the coast in HIST simulations highlight an important contribution of anthropogenic climate change to changes observed in SLP for the coastal region, which in turn is able to modify alongshore winds. The assumption of linear response to climate change comes with issues and caveats, as widely recognized (Lindzen and Giannitsis, 2002; Rodionov, 2004), including the choice of start and end points in cyclical data, and the masking of non-linear events. Nevertheless, linear regression to quantify trends over the historical record and detect responses to climate change is used frequently and remains within the state of the art of climate research (e.g., Santer et al., 2004; Boisier et al., 2016; Williams, 2017).

Trends obtained for the ERA-Interim data reveal an increase in the upwelling-favorable winds at the southern boundary of the CUS, particularly between 37°S and 42°S, while a decrease is observed at lower latitudes (30–36°S) during the austral summer. During winter months, positive trends of alongshore winds are observed between 30° and 37°S, with negative trends observed further south (Figure 3a). This seasonal spatial pattern of alongshore wind trends over south-central Chile is highly consistent with an intensification and expansion of the SPA, shifting its influence southward at the coastal region (Figure 2b). Furthermore, the spatial configuration of alongshore wind trends during the upwelling season is similar to that projected over various upwelling regions under warming scenarios. Regional projections show an increase in upwelling-favorable winds off the Chilean coast south of 35°S (Garreaud and Falvey, 2009; Goubanova et al., 2011; Belmadani et al., 2014), while GCMs show increased upwelling at the poleward boundary of the CUS. In addition, GCMs have demonstrated that the CUS off south-central Chile exhibit one of the major and more robust changes in timing, intensity and spatial distribution of upwelling-favorable winds as compared to other eastern ocean boundaries (Rykaczewski et al., 2015; Wang et al., 2015).

Several authors concur that this spatial pattern of change in alongshore winds is consistent with a poleward displacement of the SPA (Falvey and Garreaud, 2009; Belmadani et al., 2014; Rykaczewski et al., 2015) in association with a projected increase in the Southern Annular Mode and expansion of the Hadley cell under global warming (Hu and Fu, 2007; Lu et al., 2007; Seidel et al., 2008; Gillett and Fyfe, 2013). Given that the increase of upwelling-favorable winds is related to changes in SLP gradients, a notable agreement in the spatial configuration of trends in SLP gradients is observed between ERA-Interim and HIST simulation during the upwelling season (Figure S1).

However, the HIST simulation exhibits lower trends of SLP gradients, consistent with the trends of SLP presented in Figure 2. GCM simulations generally show a strengthening of upwelling-favorable winds during the austral summer at latitudes between 35°S and 40°S (Figure 3b), but the linear trend slope is generally lower than for ERA-Interim wind data. However, spatial coherence between HIST simulations and ERA-Interim data suggest that trends observed in upwelling-favorable winds over previous decades could be attributable, at least in part, to anthropogenic climate change.

A more detailed analysis of wind trends near St18 is shown in Figure S2. Accumulated Ekman transport during the upwelling months (September–February) exhibits a positive trend (p < 0.05) with some interannual variability in Ekman transport observed and partially explained by ENSO. Results show that during the negative phase of ENSO (La Niña), upwelling transport is strongest due to an increase in alongshore wind intensity; conversely, during the positive phase of ENSO (El Niño), upwelling is less intense, as shown previously by Rahn (2012).

On a decadal temporal scale, fluctuations in SST and thermocline depth have been related to the PDO (Montecinos et al., 2003; Pizarro and Montecinos, 2004). SST fluctuations in the Eastern Tropical Pacific are forced by variability in equatorial zonal winds. However, poleward of 30°S, local winds are able to reinforce the decadal variability observed in thermocline depth and SST through vertical advection associated with upwelling-favorable winds (Montecinos and Pizarro, 2005). Changes in coastal SST observed over recent decades would therefore be influenced by both PDO and the intensification of upwelling-favorable winds. In reality, the negative trend of the PDO index only partly explains the cooling trends observed remotely in SST (Falvey and Garreaud, 2009). Results show that the seasonal spatial pattern of SST trends appears to respond to wind trend patterns, particularly during the upwelling season, when the coastal ocean is experiencing clear cooling in the poleward boundary of the CUS (Figure 3c). The simulated impacts of climate change along the coast of central-south Chile include a decrease in SST of 2–3°C during the austral summer as the result of changes in alongshore winds (Aiken et al., 2011).

MLD variability is driven through a balance of processes that enhance and break down mixing and stratification (Turney and Banerjee, 2008). Momentum fluxes induced by increased winds at the air-sea interface generally result in shear and convective instabilities at the surface and lead to mixing and reduced stratification, resulting in a homogeneous surface layer. Conversely, surface heat and salt fluxes cause stratification, which can in turn enhance or reduce the MLD (Bindoff et al., 2007). At St18, estimated MLD has deepened over the period from 2002–2016, demonstrating a negative trend of –3.91 m decade^–1 (p = 0.06) for T_crt (Figure 4a), and –2.43 m decade^–1 (p = 0.07) for D_crt (Figure 4d). These trends may be due to the increase in local winds previously described, in addition to decreasing freshwater input (Garreaud et al., 2017). Additionally, seasonal analysis indicates that in winter (June–July–August) the MLD for T_crt has a negative trend of –3.85 m decade^–1 (p = 0.58; Figure 4c), and for D_crt a trend of –5.06 m decade^–1 (p = 0.44; Figure 4f), while during the upwelling season (December–January–February), the MLD trend displays a slight rise of 1.84 m decade^–1 (p = 0.25) for T_crt (Figure 4b), with a slight deepening of –1.44 m decade^–1 (p = 0.41) for D_crt (Figure 4e).

These results indicate that the winter MLD represents the greatest contribution to general deepening; they are also consistent with results obtained from simulations using a ROMS1D model based on climatological databases corresponding to the geographical location of St18 (García-Loyola S, personal communication), where sensitivity analysis forced by a variation in momentum, heat, and freshwater fluxes shows greater variation in winter as compared to summer. The slight variation observed during the austral summer can be explained by the combination of two factors: increased incoming solar radiation that heats the sea surface and generates a marked thermocline within the top few meters; and increased winds that activate the upwelling process that pushes subsurface water to the surface and produces mixing. Both of these factors act simultaneously and inhibit MLD variation. During winter months, the thermocline disappears and surface water is more prone to vertical variability.

Monthly time series (2003–2017) and annual cycles of physical (temperature and salinity) and biogeochemical variables (chl a, nutrients NO₃^– and PO₄^3–, nitrous oxide, and dissolved oxygen) are shown for St18 in Figure S3. Table S2 shows the trends obtained for monthly anomalies (from 2003 to 2017), with estimates of summer and winter trends. Temperature presents a negative trend over the entire water column, particularly within surface waters (Table S2). Temperature at St18 shows a trend of –0.13°C decade^–1 (p = 0.54) in summer and –0.68°C decade^–1 (p = 0.16) in winter. Conversely, salinity shows a positive tendency over the entire period of study, with an overall trend of 0.10 decade^–1 (p ≤ 0.05), and a winter trend of 0.20 decade^–1 (p = 0.53). Trends toward cooler temperatures and increased salinity have been reported by Schneider et al. (2017).

MLD shoaling, primarily as a result of alongshore wind patterns, may partially explain observed temperature and salinity trends. In spring–summer, the water column is pushed upwards as a result of the increased coastal upwelling, which leads to an increase in positive vertical water speed, as denser water from the subsurface is lifted to the surface. This subsurface water, mainly from the Equatorial Subsurface Water, is denser than the surface water, which is comprised principally of sub-Antarctic waters (Sobarzo et al., 2007). During spring and summer, the increase in offshore Ekman transport and Ekman pumping therefore results in denser (more saline and cooler) surface water. Additionally, the effect of reductions in precipitation and freshwater discharge from the Itata River on the winter MLD needs to be considered. The decrease of fresh water discharge has been linked to the recent megadrought affecting Chile (Garreaud et al., 2017). The freshwater layer, typically concentrated within the first 10 m, is formed during the rainy season (winter) and acts as a barrier to mixing, given the greater buoyancy of freshwater as compared to subjacent, saltier water (MacIntyre et al., 2010). The estimated trend toward increasing salinity within the top 15 m of the water column is consistent with a reduction in continental runoff and a shortening of the rainy season, resulting in reduced freshwater mixing with seawater and overall greater salinity. The winter MLD also tends to be deeper due to increased wind gusts or storm winds that drive the mixing of surface and subsurface layers. Poleward winds prevail in winter months and result in more frequent downwelling processes. It is possible that the surface cooling and denser water generated at St18 is related to a slight reduction of the MLD, which would transport a greater volume of cooler water mass towards the surface. In this way, changes in upwelling contribute to lower temperatures of surface waters.

CUS fluctuate over time and space. Several studies focused on satellite imagery for SST and chl a (Letelier et al., 2009; Aravena et al., 2014) indicate that local conditions, particularly coastline characteristics and topography, can influence upwelling (Sobarzo et al., 2016). In the short term, seasonal cycles provoke greater variability as compared with, for example, low frequency processes such as ENSO (Corredor-Acosta et al., 2015; Testa et al., 2018). Estimates of chl a trends (2003–2016) indicate some spatial heterogeneity along the coast (Figure S4). However, a positive slope is observed between 36°S and 40°S, a latitudinal section that also shows a cooling trend (Aravena et al., 2014; Schneider et al., 2017), in contrast with the 32–36°S latitudinal band, which displays a negative trend.

Ongoing monthly anomaly time-series data from St18 for surface PAR, chl a, NO₃^– and PO₄^3–, dissolved oxygen, and nitrous oxide, since 2002 are illustrated in Figure 5. Average monthly PAR and phytoplankton concentrations (chl a), derived from satellite data, exhibit positive trends: for PAR, 0.40 E m² d^–1 decade^–1, p = 0.33; for chl a, 0.37 mg m^–3 decade^–1, p = 0.38. Positive trends are also observed for dissolved oxygen and nitrous oxide, whereas nutrient concentrations decrease over time. Despite these decreases, we observed positive trends (with no statistical significance) in nitrate and phosphate fluxes toward surface waters (0.6 and 1.5 kmol d^–1 decade^–1; Figure S5), mainly explained by increases in Ekman transport (Figure S2). Table S2 displays temporal trends and statistics for these variables.

Figure 6 shows the relationship between surface satellite chl a and PAR, and between integrated chl a and nitrate flux, for the upwelling period (September–February); in both cases, the relationships are positively related by linear analysis (r² = 0.65, p < 0.05 [Fisher’s p-value] and r² = 0.41, p < 0.05, respectively). These correlations suggest greater incoming short wave radiation resulting from decreased cloud coverage produced by intense upwelling-favorable winds (Garreaud and Muñoz, 2005), along with enhanced nutrient flux toward the surface layer, which might stimulate phytoplankton growth (Wang et al., 2015).

Additionally, observed surface NO₃^– concentrations increased in summer months and decreased in winter months (Table S2). Sources of nutrient input to surface coastal waters include upwelling and continental drainage, both of which stimulate the growth of phytoplankton. With respect to freshwater discharge, an extensive megadrought has been reported in Central Chile since 2010, and has caused water scarcity, vegetation deterioration, increased forest fires, and decreases in nutrient export along the Chilean coast (Garreaud et al., 2017). Chilean rivers flowing out of south-central Chile have undergone an abrupt decrease in discharge and resulting decreases in nutrient inputs, coastal primary production, and river plumes (Masotti I, personal communication). Prior to the megadrought, a strong influence on river discharge and nutrient export to the coast over Central Chile was attributed to ENSO (Thiel et al., 2007).

These changes in nutrient flux are consistent with observations made during recent years of intensification in coastal eutrophication and land runoff (Rabalais et al., 2009), in some cases as the result of upwelling pumping (Pennington et al., 2006). When in situ measurements (integrated from surface to a 15 m depth) are considered, increased chl a is only observed during spring, revealing the ecological response to the complex changes in the phytoplankton community as a result of stratification or stability (Dave and Lozier, 2013), as is the case within the MLD (Figure 4). These in situ measurements indicate high synoptic variability that is not captured by monthly sampling.

Dissolved oxygen in surface water was observed to decrease in the spring and summer but increase in the winter (Table S2); greater advection of oxygen-poor water from the subsurface (mainly Equatorial subsurface water poor in oxygen, as described by Farías et al., 2009) to the surface is expected to be counteracted by increased oxygenation via vertical mixing and photosynthesis. For this reason, despite the fact that observed effects are consistent with those expected given seasonal trends (Table S2), understanding the net effect on trends in dissolved oxygen concentrations is difficult.

A positive trend in N₂O concentration in the mixed layer was also observed (Table S2). This increase might be driven by biogeochemical or physical processes or both. The main biogeochemical processes involved in N₂O production are microbial nitrification and denitrification. The first process operates under a wide range of oxygen concentrations, while the second takes place in suboxic and anoxic conditions (Lam et al., 2009). Although nitrification is considered to be strictly aerobic (it requires oxygen for the oxidation of ammonium and nitrite), it can persist at nanomolar oxygen levels as low as ~0.01% air saturation (Füssel et al., 2012, Bristow et al., 2016). Galan et al. (2017) recently demonstrated that most of the N₂O produced in the oxycline off central Chile comes from aerobic ammonium oxidation, while in bottom waters dissimilatory nitrite reduction becomes the (seasonally) dominant process. On the one hand, significant correlations of chl a levels with N₂O hotspots in the thermocline (oxycline), as observed in a previous study at St18, suggest that the microbial activity fueled by the increased availability of organic matter stimulates N₂O production in our study area (Farías et al., 2015; Galán et al., 2014). In a scenario of greater vertical advection of nutrients with concomitant increase in chl a, as in the present study, enhanced primary production could thus explain the stimulation of N₂O production. On the other hand, higher levels of N₂O produced in the subsurface layers could be advected and degassed quickly due to intensified wind-driven upwelling (Figure 3), leading to an increase in air-sea N₂O flux towards the atmosphere. This scenario is expected in regions with more turbulent wind regimes, where gases are rapidly transferred toward or away from the ocean surface.

It is important to note that last decades trends in sea level pressure and wind during the last decades could not be compared directly to those observed at St18, as the two data sets address different periods of time. Here, we have attributed changes to natural variability, as associated with the PDO, and to climate change using reanalysis and GCM, without any attribution using St18 data. Nevertheless, observed changes in the water column at St18 are consistent with an increase in alongshore winds, solar radiation and upwelling processes, providing insights into the effect of anthropogenic forcing on coastal upwelling. An important lesson can also be learned from the fact that the spatial pattern of alongshore wind trends during the last decades is coherent with that projected under global warming.

Under future scenarios, increases in solar radiation and wind-driven coastal upwelling of cold water might have negative repercussions on the phytoplankton community, as a stronger thermal gradient between the surface and subsurface layers may inhibit nutrient input from deeper waters into the photic zone (Boyce et al., 2010; Wang et al., 2015). Furthermore, physical constraints such as wind mixing may have favorable or adverse effects on fisheries, for example, through the mechanisms described by Cury and Roy (1989) and Bakun (1990) that lead to an “optimal environmental window” for recruitment. Decadal behavior of environmental variables may be highly variable and respond to specific thresholds of wind and/or stratification (associated with precipitation). However, whether synergic or antagonistic responses result from the multiple feedback processes between the atmosphere and the ocean surface is unclear.

Over previous decades, an increase in SLP over the Southeast Pacific has been observed, with consequences for upwelling-favorable winds, coastal SST, solar radiation, nutrients, phytoplankton biomass, and dissolved oxygen concentrations. Although natural interdecadal variability associated with the PDO likely plays a key role in changes in surface, atmospheric and oceanographic conditions, it is possible that anthropogenic climate change also contributes to or amplifies these shifts. The spatial configurations of these trends are consistent with those projected by climate change scenarios for the 21^st Century. This research demonstrates that SLP trends observed over previous decades for the coastal region of south-central Chile are not primarily forced by the PDO but more likely represent the result of anthropogenic climate change. If changes in SLP on the coast and resulting impacts on alongshore winds are attributable, at least in part, to climate change, changes observed for several related oceanographic conditions should also be at least partially attributable to climate change. As the majority of changes in surface water are attributed to changes in stratification, MLD shoaling observed over the past decade could thus be related to climate change and marked by an increase in upwelling-favorable winds and reduced precipitation and continental runoff.

The evidence summarized in this work illustrates the rapid rate of change and the complexity of the biophysical system defined by the coastal upwelling system along central and southern Chile. Observed trends are influenced by anthropogenic climate change, and the chain of impacts affects coastal ecosystems that sustain key livelihoods and, in turn, human well-being. Changing coastal upwelling systems and their consequences therefore represent a clear regional manifestation of the Anthropocene. Effectively addressing this issue will require a systems analysis approach to better integrate policy and science. To this end, it is vital to increase the level of certainty attached to the attribution of causality to past and future change by broadening the observational basis and developing long-term monitoring of the coastal upwelling system of central and southern Chile.

Gridded products and Global Climate Model simulations are openly available in several repositories. In-situ oceanographic observations at station 18 are property of the University of Conception (Chile), to get these data write to the corresponding author.

The authors acknowledge the supply of measurements of sea surface temperature provided by the Servicio Hidrográfico y Oceanográfico de la Armada de Chile (SHOA). We also thank Laura Gallardo and Rene Garreaud for their constructive comments, which helped to improve this study. We also thank G. Garcia for assistance during laboratory analyses, and the crew of the R/V Kay Kay. This study was carried out based on the COPAS (www.copas.cl) until 2012 year.

This special issue represents a contribution to the Center of Excellence FONDAP 1511009. This study received financial support from CONICYT; LF was supported by FONDECYT N°1161138, GT was supported by CONICYT-PFCHA/Doctorado Nacional/2017-21170561 and CA was supported by PAI grant N°79150062, FONDECYT N°11171163 and COSTAR-UV.

The authors have no competing interests to declare.

• Contributed to conception and design: CA, SG, GT, LF

• Contributed to acquisition of data: CA, SG, GT, DS, LF

• Contributed to analysis and interpretation of data: CA, SG, GT, DS, LF

• Drafted and/or revised the article: CA, SG, GT, LF

• Approved the submitted version for publication: CA, SG, GT, DS, LF