Among other available reanalysis products, the ERA5 product is used mainly due to its spatial resolution (0.25°×0.25°) and extended temporal coverage (1950–2020; Hersbach et al., 2020). The ERA5 surface level wind field and sea surface temperature (SST) products for the period of 1950–2019 are used. The data sets consist of two major reanalysis periods. Over 1950–1978, it was reconstructed from the assimilation of early satellite-derived data and reprocessed conventional observations (Bell et al., 2020). In contrast, vast amounts of historical observations are used after 1979, which are available a few days behind real time (Hersbach et al., 2018). ERA5 hourly products with global coverage are publicly accessible on a 0.25° × 0.25° grid. In this study, the daily mean values of wind components were computed from the hourly data.

The time series of the ENSO index (NINO 3.4) was taken from the climate prediction center of the National Oceanic and Atmospheric Administration ¹ (NOAA), which was obtained from the 3-month running mean of monthly SST. The IPO index was also taken from the NOAA Physical Sciences Laboratory (PSL) ² .

The daily τ is calculated using the bulk formula:

where ρ is the air density (1.2 kg.m⁻³), C D N is the neutral drag coefficient (Trenberth et al., 1990), and U₁₀ is the wind velocity at 10 m height above the sea surface.

Considering the impacts of wind speed on the drag coefficient is sometimes essential to accurately estimate momentum flux into the ocean over the coastal area (Brodeau et al., 2017; Bonino et al., 2022). Since the momentum flux associated with the wind is of central importance for the upwelling process in EBUSs, we used the approach of Trenberth et al. (1990), which is based on the drag coefficient varying by wind speed. This method is widely used by marine researchers.

Theory suggests that the coastal upwelling is bounded to a coastal zone in the order of the first baroclinic Rossby radius of deformation (R_i). We consider the wind over 100 km from the coastline, which is within the range of R_i across the PUS. The temporal evolution of the τ at the 100 km coastal band is calculated by zonal averaging at each latitude.

The monthly, seasonal, and annual mean of the τ and the corresponding anomaly of monthly and annual τ with the reference period of 1981–2010 are calculated. The climatological cycle of τ is calculated as well.

The anomaly of mean monthly SST is calculated using ERA5 products for 1950–2019. To estimate the wind-stress-curl-driven upwelling, we computed the curl of τ divided by the Coriolis parameter (f) as follows:

where τ_x , τ_y , f, Ω, and φ are zonal and meridional wind stresses, the Coriolis parameter, the angular speed of the earth, and the latitude, respectively. The climatological cycle for 1950–1979 and 1990–2019 is calculated. These two periods are compared to investigate the long-term trend of wind stress.

To understand the SST and surface wind field co-variability, we implement joint Empirical Orthogonal Function (EOF) analysis using ERA5 wind and SST from NOAA satellite OIv2 measurements (Kutzbach, 1967; Reynolds et al., 2007). Before conducting the EOF analysis, the climatological annual cycle was subtracted from each variable. Then, given the fact that wind and SST have different physical units and scales, we normalized them with the areal average of the standard deviation of monthly anomalies over the entire domain (85°W–70°W, 20°S–2°S). In this way, the EOF modes and their principal components (PCs) can be interpreted as the co-variability of the normalized surface wind field and SST in space and time, respectively.

The long-term mean and variability of wind stress (τ) are shown in Figure 2 . The variability is expressed as the standard deviation of daily mean wind stress. A distinct maximum τ of approximately 0.08 N/m² is observed off Pisco, which is reduced seaward. The strong τ off Pisco–Marcona is attributed to the local effect of coastal geometry and headland. The τ shows a high variation in the Pisco–Marcona area with a standard deviation of approximately 0.03 N/m². The variability is reduced across the northern and southern coasts.

The alongshore component of τ shows a large variation in seaward (i.e., zonal) direction across different cells ( Figure 2B ). Due to the headland effect of the coastline, there is a strong wind close to the coastline of the Pisco and Marcona cells, with a steep reduction toward the land. The wind off Paita shares several similarities with that off Pisco but with a smaller magnitude. This may come down to the similar orographic and headland pattern in these two upwelling cells. Unlike the central and northern portions of the HUS, τ profile shows a stepwise enhancement with the coastal distance in the sector south of Marcona.

In austral summer (January–March) concurrent with the southernmost position of the intertropical convergence zone, τ experiences its seasonal minimum in the entire PUS with values of approximately 0.02–0.05 Nm⁻² in the coastal band. The standard deviation is almost uniform, approximately 0.01–0.02 Nm^-2 in the coastal area and enhanced seaward ( Figure 1 ). The τ is more intensified in the austral autumn (April–June) and reaches the seasonal maximum in winter (July–September). Except in the southernmost portion of the HUS, the intensity of τ exceeds 0.04 Nm⁻² over the entire region in winter. It reaches approximately 0.1 Nm⁻² in the Pisco and Marcona cells. A positive seaward gradient in the τ is evident near the coast, which is more pronounced at the Pisco and Marcona. The offshore gradient becomes steeper in the Austral winter. It is important to highlight the fact that the presence of distinct maximum τ off the Pisco and Marcona cells persists throughout the year.

The annual cycle of the alongshore τ over the coastal zones with a width of R_i shows a high seasonal variation with a seasonal maximum in Austral winter (July–September) and a seasonal minimum in Austral summer (January–March). The seasonal pattern shows a non-uniform geographical distribution with a strong seasonal peak off Pisco–Marcona (13°–16°S) diminishing north- and southward ( Figure 1 ).

As mentioned in previous studies (Chavez et al., 2003; Taylor et al., 2008; Dewitte et al., 2012; Espinoza-Morriberon et al., 2017; Chamorro et al., 2018), the Pacific El Niño and La Niña events exert profound impacts on the wind field in the HUS. The mean τ during El Niño and La Niña events in 1950–2019 is calculated using monthly mean and plotted in Figure 3 . The criterion for selecting El Niño and La Niña events is the Nino3.4 SST anomalies. When the SST anomaly is higher (cooler) than >0.5°C and persistent for longer than 5 consecutive months, it is considered El Niño (La Niña). During the El Niño, the magnitude of τ shows higher values concerning the long-term mean, particularly along the coastal area, with a maximum of ~0.09 Nm⁻² in the Pisco upwelling cell ( Figures 3A, C ). On the contrary, the La Niña events are associated with reduced τ mainly in the central and northern parts ( Figures 3B, D ). The higher τ during the El Niño does not mean stronger upwelling and productivity. The deeper thermocline and thicker warmer water cause weaker upwelling and less nutrients and consequently lower primary production. On the contrary, during La Niña, the thermocline is shallower concerning normal conditions, which causes more nutrient supply to surface waters and higher primary production.

The time evolution of yearly mean alongshore τ anomalies over the 100-km coastal band is displayed in Figure 4 . It shows fluctuations on interannual to decadal time scales. A year-to-year oscillation with more signatures in the central and southern sectors of the PUS is obvious, which is likely associated with the ENSO events. Robust positive (equatorward) anomalies in 1983, 1995, 1998, 2002, 2010, and 2015 coincide with pronounced El Niño events. Moreover, a transition from negative to positive anomalies was evident around 1997. A mirror-like pattern is observed around 2015. However, the center of large positive and negative anomalies during these two events does not appear at the same location. It might be connected to the asymmetric behavior of the ENSO and the IPO. During 1997–2015, the anomaly of τ is positive over the entire coastal domain. This situation is not only dominated in the coastal domain but also appears in the offshore area (not shown here). In addition, there is a relatively weak decadal negative anomaly period for 1972–1983, which is followed by another negative in 1984–1995 after a positive signal in 1983–1984. These systematic decadal changes resemble the transition phase of the IPO.

We applied low pass filters (1- and 5-year running average) to the time series of monthly mean alongshore τ anomalies over the coastal stripe of the Peruvian coast (5S–18S) ( Figure 5 ). The presence of fluctuations in interannual and interdecadal timescales in the time series is evident. The interannual fluctuation of the wind stress is in phase with ENSO cycles. The same is true for the long-term fluctuations of τ and the IPO cycle, but they are out of phase. The pronounced shift from negative to positive τ anomalies in 1996 and from positive to negative in 2015 is apparent in the time series, which is in good agreement with IPO phase transitions. The negative phase of the IPO, which is a basin-wide La Niña-like pattern, is observed concurrent with the positive wind stress anomalies. This is also true for the period from 1975 to 2019.

The linear correlation between the time series is also examined. The correlation coefficient of monthly τ and ENSO is approximately 0.28. Considering the entire period, the correlation coefficient of monthly τ and the IPO is not statistically significant. Therefore, the time series of IPO is divided into two periods; the first period is 1950–1978, and the second period is 1979–2015. The monthly τ anomalies show a small correlation of approximately 0.09 with the first period (1950–1978). In contrast, it is anti-correlated with the second period of IPO with a correlation coefficient of approximately −0.36. The 1-year filtered data have correlation coefficients of 0.27, 0.22, and −0.57 with ENSO and with the first and second periods of IPO, respectively. The 5-year filtered data show a high negative correlation coefficient with the second period of IPO as −0.73 and 0.32 with the first period of IPO ( Table 1 ).

To investigate the co-variation of the rotation wind stress and SST, the cross-correlation between monthly mean anomalies of SST and curl(τ)/f is examined ( Figure 6 ). Since the dominant driver of upwelling near the coastal band of about the first baroclinic Rossby radius of deformation is the alongshore wind-driven upwelling (Bordbar et al., 2021), we do not show the correlation over the coastal zone. The higher value is obtained with curl(τ)/f leading by 1 month. The curl(τ)/f is proportional to the offshore vertical transport. Since the main driver of upwelling near the coast is along-shore wind-driven upwelling, we only discuss the aforementioned correlation offshore (>100km off coast). It is shown that the curl(τ)/f, and SST is widely anti-correlated: implying that the enhanced curl(τ)/f is very likely followed by SST cooling. A higher correlation is observed off the upwelling cells.

To understand the coherent variation in surface wind and SST variability, we conducted joint wind-SST EOF analysis using the monthly SST and wind anomalies ( Figure 7 ). In this way, the co-variability of the surface wind and SST is reflected in the spatial pattern of the EOFs. Here, we show the results of our analysis for the first four EOF modes. The remaining modes appear less significant, with the percentage of explained variability typically<3.5%. The first mode ( Figure 7A ) explains more than 47% of the total variability. It features northwesterly/southeasterly anomalies in the surface wind, almost the same magnitude over the entire domain. The pronounced phase transition occurred around 1995 ( Figure 7E ), which was nearly coincident with the phase transition of the IPO. The correlation of PC1 with the IPO is approximately 0.42. However, the SST variation in this mode is small relative to that for the wind. The second and third modes seem to be related to the ENSO fluctuations and explain more than 28% of the total variability ( Figures 7B, C ). The corresponding PCs display fluctuations on interannual timescales ( Figures 7F, G ). The pronounced peaks in 1982–1983, 1987–1988, 1997–1998, and 2015–2016 in the second PC reminds the strong El Niño events. The correlation between the second PC and Nino3.4 index is approximately 0.60. The wind anomaly associated with the second mode is small compared with other modes. The warming during the El Niño-like phase is associated with the local deepening of the thermocline and elevated sea level driven by upper ocean heat content transported from the western tropical Pacific. Thus, the local wind plays rather a minor role in this mode. However, the third mode, which seems to be ENSO related, suggests a cooling associated with equatorial winds with more signature near the coast. A close connection between the spatial pattern of SST and wind is viewed in the fourth mode. However, this mode explains only 4.3% of the total variability.

To investigate the influence of El Niño and La Niña on the meridional wind stress, we regressed the monthly mean wind stress onto the monthly ENSO index. It shows higher coefficients off the upwelling cells with a maximum of 0.42 in the Pisco area ( Figure 8 ). This coefficient is particularly larger off the Pisco cell. We additionally applied the regression for the El Niño and La Niña episodes separately. The coefficients are 0.61 and 0.28 for El Niño and La Niña, respectively. The El Niño and La Niña events are decomposed based on the SST anomalies over the Nino3.4 index. The anomalies larger (smaller) than 1.0 (−1.0) are considered El Niño (La Niña) events. During El Niño, wind over the coastal zones is intensified, which is more pronounced off Pisco, Lima, Chimbote, and Paita. The weakening of the wind over the coastal zone is evident during La Niña. However, the impact of La Niña appears to be not mirror-like of El Niño. For example, wind over the central and southern portion of the PUS appears to be anti-correlated with the Nino3.4 SST during La Niña. In contrast, the wind over these sectors is marginally affected by the Nino3.4 SST during El Niño. It reminds the asymmetric features of ENSO cycles. Overall, coastal upwelling is more affected during warm (i.e., El Niño) than cold (i.e., La Niña) episodes.

The long-term trend of meridional wind stress over the coastal domain and offshore areas shows a positive trend. The climatological mean over the periods of 1950–1979 and 1990–2019 and the difference between the two periods are shown in Figure 9 . The distribution of trend ( Figure 9C ) shows positive values in almost all areas except the northernmost of the domain (above 5°S) with the maximum value in the Pisco area. The long-term wind intensification across the PUS can be largely attributed to the strengthening and southwest migration of the Pacific Subtropical Anticyclone (Ancapichún and Garcés-Vargas, 2015). Detailed information on the connection between the Pacific subtropical anticyclone and wind-driven upwelling across the PUS would be very beneficial for improving regional decadal prediction, which is beyond the current study.

In addition, the trend of the meridional wind stress over the coastal zone is calculated. First, the linear trend for the austral winter (July–September) wind stress over the coastal band (100 km) is computed. Given the phase locking of the ENSO cycle, the winter-time wind stress is slightly affected by the ENSO-related oscillations. This shows a positive trend in the entire coastal domain with higher values off the Pisco upwelling cell ( Figure 10 ). The meridional wind stress shows a long-term upward trend of approximately ~0.001 Nm⁻² per decade. Second, the impact of ENSO is subtracted from the monthly wind stress using linear regression. Again, the result shows an upward trend, which is comparable with the previous trend.