Future changes in coastal upwelling and biological production in eastern boundary upwelling systems

Domain: earth_science
Published: 2024-07-24
Editorial outcome: major_revision

Abstract
Upwelling along ocean ice eastern boundaries has attracted significant attention due to its profound effects on ocean productivity and associated biological and socioeconomic implications. However, uncertainty persists regarding the evolution of coastal upwelling with climate change, particularly its impact on future biological production. Here, using a series of state-of-the-art climate models, we identify a significant seasonal advancement and prolonged duration of upwelling in major upwelling systems. Nevertheless, the upwelling intensity (total volume of upwelled water) exhibits complex changes in the future. In the North Pacific, the upwelling is expected to attenuate, albeit with a minor magnitude. Conversely, in other basins, coastal upwelling diminishes significantly in equatorward regions but displays a slight decline or even an enhancement at higher latitudes. The climate simulations also reveal a robust connection between changes in upwelling intensity and net primary production, highlighting the crucial impact of future coastal upwelling alterations on marine ecosystems.

Introduction
The coastal upwelling along the eastern boundary of ocean basins, driven by equatorward winds parallel to the coastline, facilitates the vertical transport of nutrient-rich water to the surface1,2. This nutrient influx nourishes marine ecosystems, fostering highly productive and biodiverse “hotspots”3,4, which are particularly pronounced in the four major eastern boundary upwelling systems (EBUSs): California (CalCS), Canary (CanCS), Humboldt (HCS), and Benguela (BCS) Current Systems(Fig.1a,b). Significantly, these EBUSs cover less than 1% of the global ocean surface but provide up to 7% of the global marine primary production and 20% of the world’s capture fisheries5. Hence, accurately predicting how coastal upwelling will respond to anthropogenic climate change and its subsequent impact on biological productivity is of utmost importance.

Conventionally, wind-induced Ekman transport has been regarded as the primary mechanism driving upwelling in the EBUSs6. Consequently, Bakun hypothesized that intensified land-sea thermal contrast under future warming conditions would produce stronger pressure gradients, thereby enhancing coastal upwelling in the EBUSs by fortifying alongshore winds7. However, recent studies have underscored the impact of poleward migration of the high-pressure systems on future changes in coastal upwelling, projecting an intensified upwelling in higher latitudes and weakened upwelling in lower latitudes8. Besides, under greenhouse warming, an earlier and prolonged upwelling season was also projected in the high latitudes of EBUSs, while a different trend was identified in the low latitudes9. These studies equate changes in upwelling with changes in the upwelling-favorable winds. However, a meta-analysis revealed substantial diversity in the responses of regional coastal upwelling to anthropogenic forcing, even under universally increasing wind trends in most EBUSs10–16. In particular, future changes in stratification14,15,17 and geostrophic transport16,18 could surpass the impact of winds, playing an important role in modulating the coastal upwelling. Therefore, considerable uncertainty persists regarding the evolution of coastal upwelling with climate change in terms of the intensity and timing2,9.

In addition to the modifications in upwelling itself, our current understanding of the potential impacts of its changes on productivity in EBUSs remains limited. It is generally thought that greenhouse warming should reduce the productivity in the EBUSs by inhibiting the upward nutrient supply due to intensified stratification, which is supported by a majority of studies14,15,19–21. However, a lack of consensus on future productivity projections derived from climate models challenges the notion that stratification is the primary determinant in predicting forthcoming productivity changes in the EBUSs1. Evidence exists that the coastal upwelling contributes non negligible to the long-term changes of productivity in some EBUSs22,23. Besides, potential alterations in the nitrate content of upwelling source water could also exert a significant influence on productivity23–25. However, how the evolving dynamics of upwelling affect the nutrient supply and productivity in the major EBUSs remains uncertain.

The limited observational record length introduces considerable uncertainties when analyzing the impact of global warming on coastal upwelling using historical data1. Alternatively, state-of-the-art climate models, which have outputs of vertical velocity, offer an approach to directly evaluate the response of coastal upwelling to global warming.

In this study, we employ an unprecedented high-resolution Community Earth System Model simulation (referred to as CESM-H) (see “Climate simulations” in Methods) combined with a series of current-generation climate simulations participated in the Coupled Model Intercomparison Project Phase 6 (CMIP6) (see “Climate simulations” in Methods) to evaluate the future changes of coastal upwelling and its impact on the biological productivity in the major EBUSs under a high carbon emission scenario.

Results
Seasonal advancement and prolonged upwelling season
We begin with a broad view of the climatological coastal upwelling in the EBUSs. Consistent with the existing theoretical arguments and simulations8,18,26,27, the coastal upwelling index (UI, see “Upwelling timing and intensity” in Methods) in the major EBUSs by CESM-H demonstrates significant seasonal variability. Yet, this seasonality varies within individual EBUSs depending on latitudes (Fig.1c–f). In the poleward regions (CalCS-P, CanCS-P, HCS-P and BCS-P, see “Region selection” in Methods), upwelling displays a consistent and significant seasonal cycle, with its largest value occurring in summer (June-July-August for Northern Hemisphere and December-January-February for Southern Hemisphere) and smallest values in winter. Conversely, in the equatorward regions (CanCS-E, HCS-E, and BCS-E), the seasonal upwelling cycle exhibits contrasting patterns, with maximum values in winter and minimum values in summer. The central regions (CalCS-C, CanCS-C, HCS-C, and BCS-C), situated between the poleward and equatorward regions, exhibit strong permanent upwelling but weak seasonal variability across all EBUSs. In addition, the upwelling also displays distinct seasonality in distinct seasonality in terms of timing and duration: most regions demonstrate persistent annual upwelling, while three poleward regions (e.g. CalCS-P, CanCS-P, and HCS-P) exhibit positive upwelling from early spring to fall lasting for 250–270 days, and one equatorward region (e.g. CanCS-E) displays positive upwelling from early fall to the following spring lasting for nearly 300 days (Supplementary Fig. 1). The simulated climatology of coastal upwelling in CESM-H closely aligns with that observed in the ocean reanalysis (Supplementary Fig. 2, see “Validation datasets” in Methods), albeit exhibiting a slightly stronger magnitude, potentially attributed to the higher spatial resolution of CESM-H. This agreement instills confidence in the capacity of CESM-H to project the coastal upwelling changes under greenhouse warming.

To quantitatively assess future changes in the timing of the upwelling season, we evaluate the shifting rates of upwelling onset and termination date by utilizing the daily vertical velocity output from CESM-H (see “Upwelling timing and intensity” in Methods). Figure 2a–d shows the linear trend (see “Trend analysis” in Methods) of the onset and termination date of the upwelling season for the regions with a seasonal upwelling. There is an obvious seasonal advancement and extension of upwelling season in both Hemispheres. In particular, the earlier-shifting trend of upwelling onset is most evident in the Northern Hemisphere, with a rate of 10.5 ± 4.4 days per century in the CalCS-P, 12.2 ± 2.6 days per century in the CanCS-P, and 23.7 ± 3.2 days per century in the CanCS-E, respectively. As for the Southern Hemisphere, the upwelling onset features a more minor earlier shift in the HCS-P compared to the Northern Hemisphere counterparts, with a rate of 8.0 ± 2.5 days per century. Concurrently, a later shifting in upwelling termination is also identified in these regions, which is statistically significant in CanCS-P, CanCS-E, and HCS-P. Consequently, the combination of earlier onset and delayed termination leads to an extended duration of the upwelling season in the CanCS-P, CanCS-E, and HCS-P regions, with a rate of 8.2 ± 1.2% per century (21.2 ± 3.1 days per century), 13.6 ± 1.4% per century (38.3 ± 4.0 days per century) and 6.9 ± 1.6% per century (19.3 ± 4.5 days per century), respectively. The prolonged duration of upwelling in the CalCS-P is almost entirely attributed to the advancement of upwelling onset, resulting in a more modest trend of 4.2 ± 2.7% per century or 10.8 ± 6.9 days per century. Additionally, the regions with permanent upwelling are projected to maintain persistent annual upwelling throughout the 21st century, indicating no anticipated changes in the timing and duration of the upwelling season (Supplementary Fig. 3). Notably, the upwelling seasonal cycle exhibits cross-shore discrepancies (Supplementary Fig 4), which might be attributed its diverse dynamic mechanisms23,28,29.

However, the secular changes in onset and termination dates demonstrate a similarity pattern, which is not sensitive to the width of the upwelling band selected, providing strong support for the validity of our conclusions.

To confirm whether the seasonal advance and extension of upwelling projected by CESM-H is model-dependent, we further assess the shifting rates of upwelling seasonal timing using CMIP6 models. Notably, quantifying seasonal timing changes was not feasible using the monthly output from the CMIP6 models, particularly for temporal shifts occurring within a time frame shorter than one month. Nevertheless, we were able to assess these changes by examining variations in upwelling during the transitional month (see “Shifting definition in monthly data” in Methods). Using the CMIP6 dataset, the results consistently show a significant earlier-shifting and longer-extending trend for the upwelling season both in the Pacific and Atlantic seasonal upwelling regions (Fig. 2d–f). In the Pacific basin regions, the CalCS-P and HCS-P, a majority of models (27 out of 32 for CalCS-P; 22 out of 32 for HCS-P) project an earlier onset of upwelling, with a significant proportion (19 out of 32 for CalCS-P; 23 out of 32 for HCS-P) anticipating a prolonged duration of the upwelling season. A similar scenario is observed in the Atlantic basin, where more than 72% (23 out of 32) and 78% (25 out of 32) of models successfully captured the trend of earlier upwelling occurrence in the CanCS-P and CanCS-E, respectively. Moreover, over 87% (28 out of 32) of models simulated prolonged duration of upwelling season in both CanCS-P and CanCS-E. Therefore, the evident seasonal advancement and extended duration of upwelling are not specific to CESM-H but are qualitatively reproduced by CMIP6 models.

To elucidate the potential mechanism underlying the advancement of the upwelling season, we have decomposed the upwelling during the onset month into components associated with the wind-driven Ekman transport and geostrophic transport (see “Decomposition of upwelling” in Methods). The upwelling during the onset month across all seasonal upwelling regions exhibits a consistent and significant increasing trend, indicating an earlier occurrence of upwelling events (Supplementary Fig. 5). However, the primary contributor varies among different ocean basins. In the Pacific basin, the advancement of the upwelling season is primarily attributed to changes in Ekman transport, illustrating the dominant role of wind variations. Conversely, in the Atlantic basin, geostrophic transport surpasses Ekman transport and promotes the upwelling intensification, thereby advancing the upwelling season (Supplementary Fig. 5). The enlarged Ekman transport in the Pacific basin is governed by the intensified alongshore wind stress, resulting from heightened land-sea pressure gradients (Supplementary Fig. 6). However, the pressure gradient changes are not caused by intensified land-sea thermal contrast as suggested by Bakun7, but are predominantly influenced by variations in atmospheric pressure within the ocean interior (Supplementary Fig. 6). This suggests that the intensified alongshore wind during the spring season in the Pacific basin maybe attributed to large-scale atmospheric circulation changes, potentially originating from meridional shifts of the Intertropical Convergence Zone (ITCZ) in a warming climate30.

Complex changes in upwelling intensity
The changes in upwelling timing and duration across EBUSs do not extend to the future changes in upwelling intensity, which is defined as the total volume of upwelled water (Fig. 3, see “Upwelling timing and intensity” in Methods). The trends of upwelling intensity in major EBUSs exhibit a pronounced latitudinal-dependent pattern. Specifically, the CanCS-P, located in the poleward portion of the North Atlantic, is projected to experience an increase in upwelling intensity by 25.4 ± 2.5% per century. Similarly, the CanCS-C, situated at lower latitudes, also demonstrates an increasing trend in upwelling intensity, albeit with a comparatively weaker magnitude (9.1 ± 1.3% per century) compared to the CanCS-P. In contrast, the CanCS-E, positioned in the equatorward portion, exhibits a significant negative trend (−5.8 ± 1.3% per century) in upwelling intensity. The HCS, located in the South Pacific, exhibits a similar pattern, characterized by a significant increase in upwelling intensity in its poleward region (17.1 ± 2.1% per century) and a decreasing trend of upwelling intensity in its equatorward region (−7.0 ± 1.0% per century). As for its central region (HCS-C), the upwelling intensity exhibits a weaker declining trend, with a magnitude of −1.5 ± 1.0% per century. Similar to the HCS, a decreasing trend in upwelling intensity is observed in both the equatorward and central regions of BCS, with a more pronounced magnitude in the former (−12.4 ± 0.7% per century) compared to the latter (−7.1 ± 0.6% per century). In contrast to the increasing trend observed in poleward regions in CanCS and HCS, BCS-P exhibits a slightly decreasing trend (−1.0 ± 0.8% per century) in upwelling intensity. This can be attributed to its lower-latitude location, associated with the geographical constraints imposed by the African continent. The CalCS differs from other EBUSs in that both poleward and central regions exhibit marginally declining upwelling intensity (−2.2 ± 1.8% per century for CalCS-P and −2.3 ± 2.1% per century for CalCS-C), which are statistically insignificant at the 95% confidence level. Besides, we remark that the future changes in upwelling intensity in most regions could be reproduced within the 50-km-wide coastal upwelling band, albeit with distinct magnitudes (Supplementary Fig. 7). However, large discrepancies are observed in certain EBUS, particular in the BCS-C. Specifically, the upwelling intensity in BCS-C exhibits a pronounced declining trend within the 200-km-wide band, while displaying a subtle decrease at a distance of 50 km from the coastline. This indicates that the upwelling intensity might also exhibit a divergent distribution along the cross-shore direction in BCS-C. These divergent changes in upwelling intensity could be attributed to complex factors, including the divergent changes in along-shore wind, stratification, and eddy activities12,23,28,29,31. A comprehensive understanding of the spatial variability of the upwelling changes as well as their dynamic deserves further investigation in future studies.

To further test the robustness of the trends of upwelling intensity in the EBUSs projected by CESM-H, we compare the results from the CESM-H with those models from the CMIP6. Due to the substantial storage burden, only monthly mean vertical velocity is provided by CMIP6 models. Therefore, we approximate the upwelling intensity in both CMIP6 models and CESM-H by utilizing the monthly averaged vertical velocity. However, the impact of using monthly instead of daily data does not qualitatively affect future changes in upwelling intensity, as revealed by CESM-H (Supplementary Fig. 8), giving us confidence in its reliability for projecting future changes to upwelling intensity. The projected upwelling intensity in most CMIP6 models exhibits qualitative consistency with that of CESM-H across the major EBUSs (Fig. 3e). Particularly, over 80% (26 out of 32) of models project a reduction in upwelling intensity in the CalCS. In the North Atlantic and South Pacific, more than 75% (24 out of 32) and 78% (25 out of 32) of models successfully replicate the amplification of upwelling in their poleward region and the attenuation of upwelling in their equatorward region. Regarding the BCS in the South Atlantic, models also successfully reproduce the latitudinal pattern in the secular changes of upwelling intensity: all models indicate a significant decrease in its equatorward region, with 78% (25 out of 32) showing a moderate decreasing trend in its central region, while only 53% (17 out of 32) demonstrate a weak decreasing trend in its poleward region. Thus, the projected changes in upwelling intensity vary among different EBUSs but remain consistent across different models within individual EBUSs, providing compelling evidence for its robustness. However, there are also quantitative differences among the projected trends in CESM-H and individual CMIP6 simulations, potentially attributed to the heterogeneous capacity of climate models in accurately representing the upwelling dynamics.

The response of upwelling intensity to greenhouse warming in the EBUSs is primarily attributed to changes in vertical velocity at 50 m (Supplementary Fig. 9a). Furthermore, these changes in upwelling intensity are predominantly associated with modifications in geostrophic transport across most upwelling regions (Supplementary Fig. 10), aligning with recent studies16. Specifically, in certain EBUSs such as the HCS-C, BCS-C, and BCS-P, despite projected increases in Ekman transport in the future, the overall upwelling intensity shows a decreasing trend. This reversal is caused by the strong offsetting effect of geostrophic transport, underscoring its significant influence. However, as latitudes increase, the dominant influence of geostrophic transport gradually diminishes, yielding to the prominence of wind-induced processes. The intensified Ekman transport at the poleward regions might be attributed to the poleward migration of atmospheric systems due to greenhouse warming8. However, future changes in geostrophic transport are influenced by various factors, including asymmetrical ocean warming, sea ice melting, etc.32,33. Consequently, distinct patterns of geostrophic transport are anticipated to manifest in individual EBUSs, warranting further investigation in future studies.

Methods
Climate simulations
The CESM-H simulation utilized in this study is based on version 1.3 of CESM. The model features a 0.25° resolution and 30 vertical levels for the atmosphere component and a 0.1° horizontal resolution and 62 vertical levels for the ocean component. For the oceanic component, the vertical grid space increases from 10 m near the surface to 250 m near the bottom. In particular, the thickness of each layer is identically set as 10 m in the upper 150 m, below which the thickness gradually increases with depth. The CESM-H includes a 500-year-long pre-industrial control (PI-CTRL) simulation and a 250-year-long historical and future transient simulation (HF-TNST). The HF-TNST simulation branches off from the 250th year of the PI-CTRL and incorporates historical forcings from 1850 to 2005 and representative concentration pathway 8.5 (RCP8.5) forcings from 2006 to 2100. For further details on the CESM-H simulation, please see the overview paper by Chang et al.48. Daily averaged outputs of vertical velocity at 50 meters are saved during 1920–2100.

Moreover, 32 climate models from CMIP6, which output vertical velocity, are employed in this study (Supplementary Table 1). The simulations incorporated historical forcings from 1850 to 2014 and utilized socioeconomic pathway 5-8.5 (SSP5-8.5, aligned with RCP8.5 adopted by CESM-H but accounting for socioeconomic effects) forcings from 2015 to 210049. Considerable consistency was observed between SSP5-8.5 and RCP85 scenarios50, enhancing confidence in comparing the results from CESM-H and CMIP6. In addition, 14 CMIP6 simulations with a biogeochemistry component were utilized in the analysis of primary productivity in the EBUSs (Supplementary Table 1).

Validation datasets
The Estimating the Circulation and Climate of the Ocean, Phase II (ECCO2) reanalysis product51 was utilized for the verification of the CESM-H. ECCO2 provides results stored as a 3-day mean field on a 0.25° horizontal grid with 50 vertical levels, covering the period from 1992 to the present. To facilitate comparison with CESM-H, the 3-day mean vertical velocity at 50 meters spanning from 1992 to 2022 was linearly interpolated to a daily resolution.

The observed vertically integrated primary productivity dataset is an experimental dataset harvested from the Pacific Fisheries Environmental Laboratory. The daily file is derived from satellite-derived measurements of chlorophyll-a, incident visible surface irradiance, and sea surface temperature52. The monthly outputs, comprising a composite of the one-day files, have been available from January 2003 to the present, with a spatial resolution of approximately 0.05°.

Region selection
In terms of seasonality, each EBUS is divided into three regions: the poleward region (-P), the central region (-C), and the equatorward region (-E). The poleward regions, located at high latitudes in each EBUS, experience a pronounced seasonal upwelling, with stronger upwelling predominantly occurring during summer and weaker upwelling or downwelling observed in winter. The equatorward regions, situated at low latitudes in each EBUS, also exhibit a significant seasonal upwelling. However, the upwelling reaches its peak during summer and weakens or even shifts to downwelling during winter, in contrast to the poleward regions. The central regions, positioned between the poleward and equatorward regions, exhibit a strong permanent upwelling but display limited seasonal variability across all EBUSs.

According to the aforementioned categorization principle, the CanCS is divided into three regions: the poleward region (CanCS-P), situated along the coast of the Iberian Peninsula (36°N–43°N), the central region (CanCS-C), encompassing latitudes ranging from 19°N-36°N, and the equatorward region (CanCS-E), traditionally referred to as the Mauritania-Senegalese upwelling region, spanning from 12°N to 19°N. Similarly, the HCS located in the South Pacific is also divided into three regions: HCS-E (5°S–15°S), HCS-C (15°S–36°S), and HCS-P (36°S–44°S). Likewise, the BCS situated in the South Atlantic is categorized as BCS-E (15°S-22°S), BCS-C (22°S-29°S), and BCS-P (29°S-33°S). The exception lies in the CalCS, where no winter-intensified seasonal upwelling was observed in its equatorward areas. Consequently, the CalCS region is divided solely into two regions: the central region (CalCS-C) spanning from 30°N to 40°N and the poleward region (CalCS-P) ranging from 40°N-48°N.

Upwelling timing and intensity
The upwelling timing and intensity were defined based on the upwelling index (UI), computed as the volume of upwelled water at 50 m within ~200 km (20 model grids) from the coast53.

UI(t) = \iint w(x,y,t)dxdy (1)

where w is the vertical velocity. It should be noted that the regions situated ~200 km from the coast encompass the coastal upwelling and a substantial portion of the upwelling driven by offshore wind stress curl, both of which have implications for the EBUSs’ ecosystem. It is worth noting that the mixed layer varies spatially and seasonally, and thus selecting a fixed depth as representative of the mixed layer may oversimplify its complexity. However, using a spatially varying depth poses challenges in closing the transport budget within the EBUSs27, thereby complicating the analysis of the underlying dynamics governing upwelling changes. As a compromise, we opt for a fixed depth to represent the mixed layer in this study. Nevertheless, sensitivity tests show that minor depth changes do not have a substantial impact on the results in this study, leading support to the validity of the major conclusions.

The onset of the upwelling season is the date on which the UI turns negative to positive; similarly, the termination of the upwelling season is when the UI turns positive to negative. The duration of the upwelling season can be calculated by determining the total number of days between the onset and termination dates. However, since the transition of UI could be influenced by high-frequency fluctuations54, we have applied a 1-year 4th-order Savitzky-Golay filter55 to remove the interseason signals of UI before the detection of onset and termination dates. The upwelling intensity is defined as the total volume of upwelled water during the upwelling season:

Upwelling intensity = \int_{onset}^{termination} UI(t)dt (2)

Both alterations in the duration of upwelling season and vertical velocity at a depth of 50m have the potential to modify the the upwelling intensity. To assess the impact of velocity changes on upwelling intensity changes, we recalculated the velocity-induced upwelling intensity by keeping the onset and termination dates fixed at their climatological values (1920-1949). Consequently, the duration-induced upwelling intensity is derived as the differences between the total and velocity-induced upwelling intensity.

Shifting definition in monthly data
The upwelling timing derived from monthly outputs is similar to the daily CESM-H method, albeit with an additional condition. To quantify the changes in monthly models’ outputs over time, an epoch difference was calculated from two 30-year periods, spanning 1920-1949 and 2071-2100. Changes exceeding a month would be reflected by the models directly. However, the coarse temporal resolution of models could induce a constant timing of upwelling under greenhouse warming, as it is unable to capture changes within a few days. Hence, following Brady et al.11., a positive linear trend of UI at the transition month is regarded as an advancement in the onset or delay in the termination. Similarly, it can be inferred that a negative trend of UI during the transition month reflects a later onset or an earlier termination. The changes in duration are determined by the linear trend of the sum of UI at the onset month and termination month, with a positive depicting an extending duration.

Trend analysis
The linear trends and their standard errors were obtained by regressing the time series of each upwelling metric against time using the Cochrane-Orcutt (C-O) method56. This method is employed to address the autocorrelation within the time series. The model for time series of a given quantity \theta(t) is expressed as \theta(t) = \beta_0 + \beta_1t + \epsilon(t), where the intercept (\beta_0) and linear trend (\beta_1) can be easily computed based on the ordinary least-squares method. Here, \epsilon(t) represents a stationary stochastic process with a zero mean. The standard error of the linear trend was computed using the C-O method, assuming that the autocorrelation in \epsilon(t) can be approximated as a first-order autoregressive process.

Decomposition of upwelling
The UI contributed by the wind-induced Ekman transport (UI_e) can be derived based on horizontal mass-flux divergence induced by Ekman transport27:

UI_e = \frac{1}{\rho_0 f} \tau \times k (3)

where \rho_0 is the reference seawater density, f is the Coriolis parameter varying with latitude, \tau = (\tau_x, \tau_y) is the surface wind stress, and k is the unit vector in vertical direction.

Similarly, the upwelling intensity contributed by geostrophic transport (UI_g) can be computed as the horizontal mass-flux divergence caused by geostrophic transport16:

UI_g = \int_{50m}^{0} \nabla \cdot u_g dz (4)

u_g = \frac{1}{\rho_0 f} \nabla p \times k

where u_g is the geostrophic velocity, p is the seawater pressure computed according to the hydrostatic approximation p = \rho_0 g \eta + \int_0^D \rho g dz, with g the gravity acceleration, \eta the sea surface height, and D the seawater density.

Computation of VNT
Nitrate input is an important component in nutrient supply as most of the phytoplankton growth was limited by the availability of inorganic nitrogen57. Thus, the vertical nutrient flux at 50 meters was calculated by multiplying the w by the mole concentration of nitrate (N) at the same depth58. Afterward, the VNT, similar to the upwelling intensity definition, is defined as the integral vertical nutrient transport from the onset to the termination of the upwelling season.

VNT = \int_{onset}^{termination} w(x,y,t) \cdot N(x,y,t) dxdydt (5)

To identify the effect of changes in upwelling and nutrients on VNT, we further decompose the VNT into components associated with the effects of upwelling (\Delta w), effects of nutrient (\Delta N), and their interaction effects:

\Delta w = \int_{onset}^{termination} w' \cdot \bar{N} dxdydt

\Delta N = \int_{onset}^{termination} \bar{w} \cdot N' dxdydt (6)

Interaction = \int_{onset}^{termination} w' \cdot N' dxdydt

with the overbar denoting the climatological mean value and the prime is defined as perturbations from the climatological mean value.

NPP during the upwelling season
The depth-integral primary organic carbon production over the full water column was regarded as the NPP in the ocean. Given that the euphotic layer depth is primarily shallower than 50 meters in the EBUSs59, the full-depth-integral NPP could serve as a representative indicator for the upper-ocean NPP. Subsequently, the NPP during the upwelling season is calculated as the summation of NPP starting from the onset to the termination.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This is an interesting article related to climate change and Eastern Boundary Upwelling Systems (EBUSs). Throughout the study, the authors engage with hypotheses that are currently central to the EBUS community's research. One of the major challenges is linked to the coarse resolution of global climate models (both CMIP5 and CMIP6). In this context, the authors present a high-resolution global model with notable capabilities for evaluating upwellings. Additionally, they introduce an innovative methodology for accurately representing upwelling seasons. Finally, a key uncertainty regarding climate change in EBUSs is the fate of nutrients under climate change conditions. In this section, the authors address this uncertainty. The article is well-written and well-structured.\nHowever, in my opinion, the results obtained in this study should undergo a thorough review before being published. My perspective is primarily supported by the high resolution demonstrated by the global model proposed by the authors. While one of the main uncertainties of EBUSs could be addressed with a high-resolution model, the authors negate this advantage by using averages at 200 km offshore, failing to adequately reflect the importance of regions on the shelf in EBUSs. Additionally, the authors do not accurately represent vertical transport in upwelling, as they use arbitrary depths for the four EBUSs with their respective latitudes and seasons. Alongside this, the authors present initial hypotheses about climate change and EBUSs, display relevant results (which they do not discuss) regarding the increase or decrease in upwelling but do not discuss or investigate which of the initial hypotheses align with their findings. Moreover, they introduce in the discussion the relationship of nutrients with stratification, when no results of oceanic stratification have been presented. Finally, the nutrient results are solely studied based on the CMIP6 global models, diminishing the significant advantage offered by the high-resolution model.\nTherefore, I suggest that the authors undertake a comprehensive analysis of the comments I have provided below. Thus, the article should be resubmitted after a major revision.\nDetailed comments:\nL 31. Add that the transport is vertical towards the surface.\nL 36 In Fig. 1 and subsequently in the rest of the sample analysis, the four EBUSs are shown. However, in the C-ICS, the authors present a domain cut off at approximately 26ºN, where upwelling is most intense throughout the year (21ºN -26ºN; Vázquez et al. 2022), and the southernmost region (Mauritania-Senegalese) does not appear, which is one of the most intense in the C-ICS (Sylla et al. 2019). Additionally, the P-CCS region presents annual permanent upwelling from 5ºS to 40ºS, where it becomes seasonal as shown by the authors. In this case, the authors decide not to conduct the analysis from 17ºS to 40ºS (Kämpf and Chapman, 2016). Therefore, I propose that the authors adjust the domains well to have a more exhaustive examination of their results and subsequent conclusions.\nVázquez, R., Parras-Berrocal, I., Cabos, W., Sein, DV., Mañanes, R. and Izquierdo, A. (2022). Assessment of the canary current upwelling system in a regionally coupled climate model. Clim. Dyn. 58, 69–85. https://doi.org/10.1007/s00382-021-05890-x.\nSylla, A., Mignot, J., Capet, X., Gaye, A.T., 2019. Weakening of the Senegalo–Mauritanian upwelling system under climate change. Clim. Dyn. 53 (7–8), 4447–4473. https://doi.org/10.1007/s00382-01904797-y.\nKämpf, J., Chapman, P. (2016). The Peruvian-Chilean Coastal Upwelling System. In: Upwelling Systems of the World. Springer, Cham. https://doi.org/10.1007/978-3-319-42524-5_5\nL41-42 The reference to the lack of observational data is from 1973. Observational systems have advanced since then. Please, look for a more recent reference.\nL43 Since the authors mention Bakun's hypothesis, please mention to the second hypothesis proposed in EBUSs’s environment by Rykaczewki et al. (2015).\nRykaczewski, R.R., Dunne, J.P., Sydeman, W.J., Garcia-Reyes, M., Black, B.A., Bograd, S. J., 2015. Poleward displacement of coastal upwelling-favorable winds in the ocean’s eastern boundary currents through the 21st century. Geophys. Res. Lett. 42, 6424–6431. https://doi.org/10.1002/2015GL064694.\nL48-50 Neither of the two provided references is related to the author's proposed statement. I suggest to the authors to search for information and modify the citations. Here's a suggestion in this regard:\nWang, D.W., Gouhier, T.C., Menge, B.A., Ganguly, A.R., 2015. Intensification and spatial homogenization of coastal upwelling under climate change. Nature 518. https://doi.org/10.1038/nature14235\nL68-69 What is the vertical resolution of the model? Please, also specify if this resolution is constant in the vertical and what is the thickness of each layer (especially in the first 200 m of the water column).\nL79-84 It would help the understanding of the manuscript if you change austral and boreal summer-winter to the respective months they refer to (e.g., JJA – June-July-August).\nL79 Please check the units for Figure 1c-f, as if the UI defined in the methodology is used, the units should be m3 s-1.\nL89 I have a question about the comparison period between the model and the reanalysis. If the model is driven under the RCP8.5 scenario, this simulation begins in 2005. Therefore, using the period 1992-2022 averages both historical model data and RCP8.5 scenario data. Wouldn't it be a better analysis to compare the period, for example, from 1992 to 2005?\nL89 Regarding the reanalysis, ECCO2 has a horizontal resolution of 0.25ºx0.25º, which is coarse for comparison with an ocean model of such high resolution as CESM-H. I propose that you perform a comparison with the GLORYS model from Copernicus, which has a resolution of 0.083ºx0.083º and may better align with the model.\nDrévillon M, Regnier C, Lellouche JM, Garric G, Bricaud C, Hernandez O (2018) CMEMS-GLO-QUID-001-030, 1.2 edn. E.U.Copernicus Marine Service Information [Online]. https:// resources. marine. coper nicus.eu/docum ents/QUID/CMEMS- GLOQUID-001- 030.pdf\nL87-91 If I am not mistaken, the model overestimates the UI compared to the reanalysis data (supplementary Fig. 1 and Fig. 1). Please add a brief sentence discussing this fact. Although the differences are within an acceptable range, it is important to note in the manuscript to encourage the modelling community to strive for improvements in the future.\nL94 I find the methodology proposed by the authors very interesting. However, I have a question about how you calculate the onset and termination of upwelling. Typically, the change from positive to negative in upwelling sign can be influenced by local wind patterns or even mesoscale processes from one day to the next. Have the authors applied any filters or analyses to calculate the average sign change over several days? Or have you simply calculated the sign change from one day to the next, and this does not vary from onset to termination?\nL95 I have a question about the choice of poleward regions. In this analysis, the Benguela upwelling does not appear as the authors have not subdivided this system. Why? The BCS has the same latitudinal extent as the California upwelling, and the seasonal latitude-dependent difference is distinguishable in Fig. 1f. For example, it could be divided at 27ºS. In this regard, please also include the BCS in the rest of the analyses, such as in Fig. 2.\nL95. Add information to the graph to clarify that the y-axis of Fig. 2a-c represents dates.\nL117. The study is conducted with the CESM-H model forced by the RCP8.5 scenario from CMIP5, but the comparison with global models is from CMIP6 (it is not mentioned in the text whether it is under the SSP5-8.5 scenario). Are there differences between the different scenarios? Add a sentence to the text or methodology discussing this.\nL133 Change Pacific to Atlantic\nL140-146 Please add the Benguela subdivision to compare whether it also increases, as in the Chile region.\nL169-173 I do not believe this study provides results to claim that it can assess local nutrient effects. By averaging over hundreds of kilometers both latitudinally and longitudinally, the ability to evaluate local effects is compromised (see Lovecchio et al. 2017, 2018; Santana-Falcón et al. 2020). Please avoid mentioning that local changes in nutrients are being assessed.\nLovecchio E, Gruber N, Münnich M, Lachkar Z (2017) On the longrange offshore transport of organic carbon from the Canary Upwelling System to the open North Atlantic. Biogeosciences 14(13):3337–3369. https:// doi. org/ 10. 3929/ ethz-b- 00019 0480\nLovecchio E, Gruber N, Münnich M (2018) Mesoscale contribution to the long-range offshore transport of organic carbon from the Canary Upwelling System to the open North Atlantic. Biogeosciences 15:5061–5091. https:// doi. org/ 10. 5194/ bg- 15- 5061- 2018\nSantana-Falcón Y, Mason E, Arístegui J (2020) Offshore transport of organic carbon by upwelling filaments in the Canary Current System. Prog Oceanogr 184(April):102322. https://doi.org/10.1016/j.pocean.2020.102322\nL173-174 I understand that the CESM-H model lacks biogeochemistry, so the authors use the global CMIP6 models to test the hypothesis. Although the results obtained previously by CMIP6 models could be compared with those of CESM-H, no evaluation of the models has been provided to see if they resolve latitudinal upwelling as in the case of CESM-H. In fact, there are different studies (Varela et al. 2022; Sylla et al. 2022) showing how CMIP6 models differ from observations. Therefore, I believe this section provides useful information that complements previous results and helps take a first step regarding the biogeochemical aspect of upwelling. However, a paragraph should be added discussing the deficiencies of CMIP6 models in reproducing EBUSs, primarily associated with the lack of resolution.\nVarela, R.; DeCastro, M.; Rodriguez-Diaz, L.; Dias, J.M.; Gómez-Gesteira, M. Examining the Ability of CMIP6 Models to Reproduce the Upwelling SST Imprint in the Eastern Boundary Upwelling Systems. J. Mar. Sci. Eng. 2022, 10, 1970. https://doi.org/10.3390/jmse10121970\nSylla, A., Sanchez Gomez, E., Mignot, J., and López-Parages, J.: Impact of increased resolution on the representation of the Canary upwelling system in climate models, Geosci. Model Dev., 15, 8245–8267, https://doi.org/10.5194/gmd-15-8245-2022, 2022.\nL185-186 The mentioned section in the methods does not exist.\nL216-219 Do you compare the Mixed Layer Depth (MLD) and nutrients based on averages from the areas shown in Fig 1? If so, I believe changes in MLD may be inadequately represented as it can be highly variable from one grid cell to another. Additionally, the authors do not assess stratification at any point. I suggest they calculate the Brunt-Väisäla frequency as an indicator of stratification changes (see Vázquez et al. 2023)\nVázquez, R., Parras-Berrocal, I. M., Koseki, S., Cabos, W., Sein, D. V., & Izquierdo, A. (2023). Seasonality of coastal upwelling trends in the Mauritania-Senegalese region under RCP8. 5 climate change scenario. Science of the Total Environment, 898, 166391.\nMajor comments:\nThe first question I have regarding the results in this work is related to the use of model vertical velocity outputs, as well as the use of an arbitrary depth (50 m) as the most representative of upwelling. I would like the authors to justify why they chose to use 50 m. In studies on upwelling systems, it has been agreed upon that there is significant variability both latitudinally and seasonally in the upwelling source depth. Using 50 m as a criterion oversimplifies the complexity of EBUSs. In this regard, using the MLD as the depth at each latitude and for each season would be more appropriate. However, estimating the appropriate depth remains highly complex. In fact, Jacox et al. (2018) proposes another method to assess upwelling volume that is superior to model vertical velocity outputs (Marchesiello and Estrade, 2010; Rossi et al. 2013; Oerder et al. 2015). In this context, I include the paragraph that describes this issue:\n“However, we do not use modeled vertical velocities to construct the index for several reasons: First, the calculation of upwelling using modeled vertical velocity is not straightforward. One must first define a representative depth at which vertical velocity is extracted (e.g., the MLD), which will be different for each grid cell, making it difficult to close the transport budget. For example, when the MLD differs between adjacent grid cells, water can enter/exit the surface mixed layer horizontally and that component of the transport will be missed from the upwelling/downwelling estimate.”\nTherefore, I would like to suggest to the authors that they reconsider the use of the MLD as the depth criterion and discuss the use of vertical velocity versus the wind field.\nMarchesiello, P., and P. Estrade (2010), Upwelling limitation by onshore geostrophic flow, J. Mar. Res., 68, 37–62\nRossi, V., Feng, M., Pattiaratchi, C., Roughan, M., & Waite, A. M. (2013). On the factors influencing the development of sporadic upwelling in the Leeuwin Current system. Journal of Geophysical Research: Oceans, 118(7), 3608-3621.\nOerder, V., Colas, F., Echevin, V., Codron, F., Tam, J., Belmadani, A., 2015. Peru-Chile upwelling dynamics under climate change. J. Geophys. Res. Oceans 120 (2), 1152–1172. https://doi.org/10.1002/2014JC0102, 99.\nJacox, M.G., Edwards, C.A., Hazen, E.L., Bograd, S.J., 2018. Coastal upwelling revisited: Ekman, Bakun, and improved upwelling indices for the U.S. West Coast. J. Geophys. Res. Oceans 123, 7332–7350. https://doi.org/10.1029/2018JC014187.\nThe second concern is related to the domain chosen for the zonal average. In this regard, the authors average over 200 km offshore and mention that this is most appropriate for capturing both Ekman transport-associated upwelling and upwelling associated with Ekman transport divergence. Averaging over such a wide portion of the coast may obscure the upwelling transport. Perhaps with low-resolution global models, obtaining a 200 km average is necessary. However, with a model of such high resolution as proposed in this study, the authors should consider narrowing this region to a maximum of 100 km. This would allow for a more appropriate assessment of upwelling, where mesoscale processes play a significant role. Therefore, I propose that the authors conduct the analysis by limiting the average to the first 100 km of the coast. In fact, in the methodology, the authors mention the work of Du et al. (2023), who also conducted a study with an average over the first 50 km offshore. Moreover, in Fig. 1a of the manuscript, the vertical velocity associated with upwelling is closely confined to the coast.", "reviewer_2": "This paper is a model-based study (CESM-H, CMIP6, ECCO2) aimed at understanding coastal upwelling variability in the context of the climate change. While the authors have undertaken significant efforts to comprehend the seasonal progression and upwelling intensity in the EBUSs, there are concerns regarding certain processes and model biases than prevent me for approving its publication.\nMajor concerts:\n1) Recently, it has been pointed out that the geostrophic flow may play a significant role in controlling the upwelling trends under the climate change (Jing et al., 2023). The authors should include an explanation regarding whether and why it is not being considered in this study. Could the results be the same if the geostrophic flow is considered?\n2) For the computation of the upwelling timing and intensity, the vertical velocity at 50 m depth has been considered. Could the results be different using another depth? (Note that Xiu et al. (2018) shows that, at least for CCS, there are three major drivers of the vertical velocity changes which are depth-depending: alongshore winds, wind-curl and mesoscale activity)\n3) How do the VNT and NPP appear under historical period compared to observational data (e.g., 2003-2022)? Most of the models exhibit a significant bias in the EBUSs, particularly with the biogeochemical variables. Can we have confidence in the results considering these biases?\nMinor points:\nLine 111-114: Is this valid for Chilean and Peruvian region? (as far as it is known, these two regions have different upwelling regimes, which one being permanent and the other seasonal)\nLine 235 Fig.1: why are figures 1c and 1d composed of two subfigures? (Should they be like 1f?)\nLine 515 Supp. Fig.2b and 2c: the curves in these figures represent the UI only for one of the initially selected subdivisions for each upwelling system (as shown in figure 1)?\nLine 547 Supp. Table 1: Include units on the table\nReferences:\nJing, Z., Wang, S., Wu, L. et al. Geostrophic flows control future changes of oceanic eastern boundary upwelling. Nat. Clim. Chang. 13, 148–154 (2023). https://doi.org/10.1038/s41558-022-01588-y\nXiu, P., Chai, F., Curchitser, E.N. et al. Future changes in coastal upwelling ecosystems with global warming: The case of the California Current System. Sci Rep 8, 2866 (2018). https://doi.org/10.1038/s41598-018-21247-7", "reviewer_3": "Title: Review of ‘Future changes in coastal upwelling and biological production in eastern boundary upwelling systems’\nThe paper investigates future changes in the duration and intensity of upwelling in the eastern boundary regions and its potential consequences on primary production. Utilizing simulations from CMIP6 models and high-resolution CESM, the authors showed that the duration of upwelling is projected to be longer under high-emission scenario. But changes in the upwelling intensity show region-specific features. The authors further explore the implications of these intensity changes on net primary production. While the results are interesting and seem to be scientifically sound, this paper still needs substantial work before it can be published. I summarize my concerns and comments below.\nMajor concerns:\n1) Refinement of introduction. This paper introduces a novel investigation into changes in upwelling duration, a distinctive aspect that warrants attention. However, it is crucial to incorporate relevant findings from existing literature to enrich the contextual background. For instance, Wang et al. (2015), utilizing CMIP5 models, observed an increased duration of upwelling in the high latitudes of major EBUSs, with differing trends in low latitudes. Summarizing these findings explicitly in the introduction is essential to provide readers with a comprehensive understanding of the advancements in this field. In addition, the introduction of future changes in net primary production only discussed ocean stratification and upwelling intensity. However, nutrients in the source water can also impact future changes in net primary production (Rykaczewski and Dunne, 2010). Please make the introduction more comprehensive.\n2) Refinement of Chile and Canary-Iberian upwelling systems definition. The definition of the Chile and Canary-Iberian upwelling systems requires careful consideration based on the results presented in Figure 1a. The depicted large upward vertical velocity extends beyond the designated boxes in the Chile and Canary regions, reaching substantially equatorward. Notably, the surface wind in the Chile box is perpendicular to the coast, which may not accurately represent the typical Chilean upwelling regions. Additionally, it is observed that stronger coast-parallel winds (as shown in Figure 1a) and higher net primary production (as indicated in Figure 1b) occur south of the defined Canary box.\nIn addition, if authors define upwelling regions based on vertical velocity, Figure 1a demonstrates that positive vertical velocity has very narrow zonal extents. This suggests that selecting 200 km away from the coast will include regions with negative vertical velocity. Given that upwelling timing and intensity calculations rely on box averages, the selection of these boxes has the potential to influence the results of future changes in upwelling intensity and duration. To ensure the accuracy of the study's findings, a reconsideration of the chosen boxes for defining the Chile and Canary-Iberian upwelling systems is recommended. This adjustment should aim to better capture the representative characteristics of these upwelling regions and minimize potential biases in the assessment of future changes in upwelling dynamics.\n3) Lack of mechanisms in upwelling changes. One notable limitation in the current study lies in the absence of a detailed exploration of the underlying mechanisms driving changes in the duration and intensity of upwelling. This omission weakens the novelty of the study, especially when compared to the work of Wang et al. (2015), who utilized CMIP5 models and already delved into discussions on these changes, demonstrating general consistency with the results presented in this paper. To strengthen the scientific contribution and novelty of this study, it is imperative to incorporate a thorough examination of the mechanisms influencing the observed changes in upwelling duration and intensity.\n4) Refinement of methodology concerning shifting definition in monthly data. A critical aspect requiring attention in the methodology is the method of shifting definition employed in analyzing upwelling duration based on monthly data, utilizing the epoch difference technique (Brady et al., 2017). While this method proves effective for studying the emergence of anthropogenic impacts on upwelling intensity, its applicability to upwelling duration analysis is questionable. As shown in Figure 2a-c, the onset and termination time suffer strong interannual/decadal variabilities. From 1920 to 2010, onset time varies from 01-15 to 03-01 in northern CalCS and Iberian system, and from 07-15 to 08-15 in Chilean system. Such variability suggests that changes in onset time using a 30-year epoch difference can be influenced by internal variability rather than anthropogenic forcing, especially when dealing with monthly data. The same problem arises from the changes in termination time. This implies that the comparison between CMIP6 and CESM-H, based on this methodology, may not be a fair one. To address this issue and ensure a more robust comparison, it is advisable to focus on confirming the role of surface wind in influencing changes in upwelling duration. Subsequently, employing daily wind data to estimate duration in CMIP6 models, as demonstrated by Wang et al. (2015), can provide a more accurate representation of the temporal dynamics of upwelling.\nMinor concerns:\n1) Descriptions of upwelling regions are incomplete. Firstly, the authors did not explicitly define Northern CalCS, Iberian, and Chilean regions, potentially confusing readers. Secondly, the authors should describe methods to calculate the trend errors of onset and termination time that are shown in Figure 2a-c. Thirdly, the term \"equatorward regions\" is used without a precise definition.\n2) What are the x-axis and numbers -3, -2, ….2, 3 in Figure 2d-f? Authors should describe it in the caption.\n3) The definition of upwelling intensity in the study is based on the time and spatial integral of vertical velocity at 50 m, with onset and termination times determined using monthly data. However, an important aspect requiring clarification is the attribution of changes in upwelling intensity. The study does not explicitly specify whether alterations in upwelling intensity are primarily driven by variations in upwelling duration, vertical velocity at 50 m, or a combination of both factors.\n4) The authors consistently express the agreement percentages between CMIP6 simulations and CESM-H, but it is crucial to note that not all CMIP6 simulations were utilized in the paper. To enhance precision and clarity in reporting, it is recommended to present the results in a format that specifies the number of simulations in agreement, such as \"xx out of xx simulations.\" This format provides a more accurate representation of the agreement statistics, ensuring that readers have a clear understanding of the scope and basis for the reported percentages.\n5) CESM-H models used in this study do not include biogeochemistry components, which means that nitrate in CESM-H is a passive tracer. Given this, it is not surprising that changes in vertical velocity dominate changes in vertical nitrate fluxes, as no additional sources and sinks of nitrate are considered in the future. To deepen the understanding of the biogeochemical dynamics, it is recommended to explore whether, in CMIP6 models with biogeochemistry components, changes in vertical nitrate flux are similarly dominated by changes in vertical velocity. This comparison will shed light on the role of biogeochemistry in influencing vertical nitrate flux patterns. In addition, total vertical nitrate flux calculated using monthly data does not include the component induced by submonthly variabilities (equivalent to eddy component). The eddy flux is resolved by CESM-H but parameterized in 1° or coarser CMIP6 models. Therefore, it is advisable to quantify the eddy fluxes, leveraging the advantages of CESM-H, to provide a more comprehensive understanding of the nuances in vertical nitrate fluxes and their implications.\nSpecific comments:\nL34: EBUSs generally stands for eastern boundary upwelling systems, not coastal upwelling systems.\nL41-42: This statement is not accurate. The reason for using wind to estimate upwelling is because upwelling is mainly wind-driven, not because sparsity of oceanic observations.\nL47-48: Previous studies showed that projected changes in wind and upwelling both show inconsistency across different models. But in one certain model, if it is not wind determining future changes in upwelling intensity, what are the other factors? Please clarify it here.\nFigure 1a: vertical velocity at which depth?\nL76: ‘Consistent with the existing theoretical arguments and simulations’, please add references here.\nL83: The maximum between 14°S-17°S in P-CCS and 27°S-28°S in BCS occurs in the early austral spring. Please precisely define poleward and equatorward regions for P-CCS and BCS.\nL86-87: Because duration of upwelling is the major focus of this paper, it will be more comprehensive to provide a figure including duration as a function of latitude in each upwelling system. It is difficult to follow the authors’ argument without any figures like that.\nL104: ‘all EBUSs’ means Northern CalCS, Iberian, and Chilean regions or regions defined by boxes in Figure 1a?\nL129: Authors divided EBUSs into poleward and equatorward regions, which show different changes in duration (Figure 2 and Supplementary Figure 2). ‘The consistent changes across all EBUSs …’ is confusing.\nL123-124: 75% and 73% are equivalent to how much out of how much?\nL133: ‘North Atlantic’?\nL134-135: Canary and Iberian regions are not defined.\nL138: north and south portions of CalCS are not defined.\nL140-141: As shown in Figure 1, the Benguela system extends from 33°S to 15°S. Authors cannot simply group it into equatorward systems.\nL157, 159, 161: Please add actual numbers (xx out of xx) after the percentage because you are not using all models that participated in CMIP6.\nL336: Should be ‘Brady et al.’\nReferences:\nWang, D., Gouhier, T. C., Menge, B. A., & Ganguly, A. R. (2015). Intensification and spatial homogenization of coastal upwelling under climate change. Nature, 518(7539), 390-394.\nBrady, R. X., Alexander, M. A., Lovenduski, N. S., & Rykaczewski, R. R. (2017). Emergent anthropogenic trends in California Current upwelling. Geophysical Research Letters, 44(10), 5044-5052.\nRykaczewski, R. R., & Dunne, J. P. (2010). Enhanced nutrient supply to the California Current Ecosystem with global warming and increased stratification in an earth system model. Geophysical Research Letters, 37(21)."}, "author_response": "Reply to the first reviewer\nWe are very grateful to you for your time in carefully reading our manuscript and\nproviding helpful comments that make our manuscript better. We have carefully\nconsidered each of your comments (in blue) and revised the manuscript accordingly.\nPlease find our response (in black) to your comments below.\nReviewer #1 (Remarks to the Author):\nThis is an interesting article related to climate change and Eastern Boundary Upwelling\nSystems (EBUSs). Throughout the study, the authors engage with hypotheses that are\ncurrently central to the EBUS community's research. One of the major challenges is\nlinked to the coarse resolution of global climate models (both CMIP5 and CMIP6). In\nthis context, the authors present a high-resolution global model with notable capabilities\nfor evaluating upwellings. Additionally, they introduce an innovative methodology for\naccurately representing upwelling seasons. Finally, a key uncertainty regarding climate\nchange in EBUSs is the fate of nutrients under climate change conditions. In this section,\nthe authors address this uncertainty. The article is well-written and well-structured.\nHowever, in my opinion, the results obtained in this study should undergo a thorough\nreview before being published. My perspective is primarily supported by the high\nresolution demonstrated by the global model proposed by the authors. While one of the\nmain uncertainties of EBUSs could be addressed with a high-resolution model, the\nauthors negate this advantage by using averages at 200 km offshore, failing to\nadequately reflect the importance of regions on the shelf in EBUSs. Additionally, the\nauthors do not accurately represent vertical transport in upwelling, as they use\narbitrary depths for the four EBUSs with their respective latitudes and seasons. Alongside this,\nthe authors present initial hypotheses about climate change and EBUSs, display\nrelevant results (which they do not discuss) regarding the increase or decrease in\nupwelling but do not discuss or investigate which of the initial hypotheses align with\ntheir findings. Moreover, they introduce in the discussion the relationship of nutrients\nwith stratification, when no results of oceanic stratification have been presented. Finally,\nthe nutrient results are solely studied based on the CMIP6 global models, diminishing\nthe significant advantage offered by the high-resolution model.\nTherefore, I suggest that the authors undertake a comprehensive analysis of the\ncomments I have provided below. Thus, the article should be resubmitted after a major\nrevision.\nDetailed comments:\nL 31. Add that the transport is vertical towards the surface.\nWe have revised it as “…, facilitates the vertical transport of nutrient-rich water to the\nsurface.”\nL 36 In Fig. 1 and subsequently in the rest of the sample analysis, the four EBUSs are\nshown. However, in the C-ICS, the authors present a domain cut off at approximately\n26ºN, where upwelling is most intense throughout the year (21ºN -26ºN; Vázquez et al.\n2022), and the southernmost region (Mauritania-Senegalese) does not appear, which is\none of the most intense in the C-ICS (Sylla et al. 2019). Additionally, the HCS region\npresents annual permanent upwelling from 5ºS to 40ºS, where it becomes seasonal as\nshown by the authors. In this case, the authors decide not to conduct the analysis from\n17ºS to 40ºS (Kämpf and Chapman, 2016). Therefore, I propose that the authors adjust\nthe domains well to have a more exhaustive examination of their results and subsequent\nconclusions.\nThanks for your comment. In the revised manuscript, we have separately subdivided\nthe Canary Current, Humboldt Current and Benguela Current systems into three regions\n(equatorward region, central region and poleward region) based on differences in\nupwelling seasonality. In this case, the Mauritania-Senegalese is encompassed within\nthe equatorward region of the Canary Current system and is analyzed separately.\nSimilarly, the permanent upwelling region from 15ºS to 40ºS is encompassed within\nthe central Humboldt Current system and also analyzed in the revised manuscript. The\nsubsequent results and conclusions were revised accordingly, with the elaborate\ndefinition of the regions added in line 422-443.\nL41-42 The reference to the lack of observational data is from 1973. Observational\nsystems have advanced since then. Please, look for a more recent reference.\nThanks for your comments. Following other reviewers’ comment, this sentence has\nbeen removed.\nL43 Since the authors mention Bakun's hypothesis, please mention to the second\nhypothesis proposed in EBUSs’s environment by Rykaczewki et al. (2015).\nThanks for your comments, we have added the hypothesis by Rykaczewki et al. (2015)\nto the introduction, see in line 45-48\n“However, recent studies have underscored the impact of poleward migration of the\nhigh-pressure systems on future changes in coastal upwelling, projecting an intensified\nupwelling in higher latitudes and weakened upwelling in lower latitudes (Rykaczewski\net al., 2015).”\nL48-50 Neither of the two provided references is related to the author's proposed\nstatement. I suggest to the authors to search for information and modify the citations.\nHere's a suggestion in this regard:\nWang, D.W., Gouhier, T.C., Menge, B.A., Ganguly, A.R., 2015. Intensification and\nspatial homogenization of coastal upwelling under climate change. Nature\n518. https://doi.org/10.1038/nature14235\nThanks for your suggestion, the reference has been changed to:\nGarcía-Reyes, M., Sydeman, W. J., Schoeman, D. S., Rykaczewski, R. R., Black, B. A.,\nSmit, A. J., & Bograd, S. J. (2015). Under pressure: climate change, upwelling,\nand eastern boundary upwelling ecosystems. Frontiers in Marine Science, 2.\nRetrieved from https://www.frontiersin.org/articles/10.3389/fmars.2015.00109\nWang, D., Gouhier, T. C., Menge, B. A., & Ganguly, A. R. (2015). Intensification and\nspatial homogenization of coastal upwelling under climate change. Nature,\n518(7539), 390–394. https://doi.org/10.1038/nature14235\nL68-69 What is the vertical resolution of the model? Please, also specify if this\nresolution is constant in the vertical and what is the thickness of each layer (especially\nin the first 200 m of the water column).\nThanks for your comment. For the oceanic component of CESM-H, there are 62 levels\nin the vertical with increasing grid space from 10 m near the sea surface to 250 m near\nthe bottom. In special, the thickness of each layer is identically set as 10 m in the upper\n150 m, below which the thickness is gradually increase with depth. We have added the\ndetailed information on the model in the revised manuscript, see lines 391-394.\nL79-84 It would help the understanding of the manuscript if you change austral and\nboreal summer-winter to the respective months they refer to (e.g., JJA – June-July-\nAugust).\nThanks for your suggestion. Since we have changed the expressions of the upwelling\nseasonality in the revised manuscript, thus, the usage of austral/boreal summer/winter\nhas been replaced to summer/winter. However, we have still incorporated the\ncorresponding months for the summer and winter in the Northern and Southern\nHemispheres into the revised manuscript, please see in line 91-92\nL79 Please check the units for Figure 1c-f, as if the UI defined in the methodology is\nused, the units should be m3 s-1.\nSorry for our mistake, this has been revised. See Figure 1.\nL89 I have a question about the comparison period between the model and the\nreanalysis. If the model is driven under the RCP8.5 scenario, this simulation begins in\n2005. Therefore, using the period 1992-2022 averages both historical model data and\nRCP8.5 scenario data. Wouldn't it be a better analysis to compare the period, for\nexample, from 1992 to 2005?\nAs reported in previous studies, the upwelling in the Eastern Boundary Upwelling\nSystems is influenced by natural variability and exhibits inherent multi-scale temporal\nvariability. Therefore, selecting a longer period for comparing the simulated vertical\nvelocity with observations would enable the removal of natural variability-induced\ninfluences. Traditionally, employing a minimum 30-year average is considered\neffective in mitigating the natural variability and widely used in climate studies (Ding\net al., 2021; Rykaczewski et al., 2015). We pretty agree with you that employing the\nRCP simulation may introduce uncertainties when compared to the observations.\nHowever, considering that current CO emissions estimates align closely with the\nmodel-reanalysis comparison is considered reasonable.\nWe have also performed a sensitivity analysis to evaluate the influence of selecting\ndifferent time period for comparing the simulated vertical velocity with the\nobservations. As shown in Fig. R1 and R2, the CESM-H also exhibit a high level of\nconsistency with ECCO2 during 1992-2005, albeit with some minor quantitative\nvariations.\nFigure R1. Comparison between CESM-H (a-d) and ECCO2 (e-h) during 1992 to 2005.\nFigure R2. Same as Figure R1 but for 1992 to 2022.\nDing, H., Alexander, M. A., & Jacox, M. G. (2021). Role of Geostrophic Currents in\nFuture Changes of Coastal Upwelling in the California Current System.\nGeophysical Research Letters, 48(3), e2020GL090768.\nhttps://doi.org/10.1029/2020GL090768\nPeters, G. P., Andrew, R. M., Boden, T., Canadell, J. G., Ciais, P., Le Quéré, C., et al.\n(2013). The challenge to keep global warming below 2 °C. Nature Climate Change,\n3(1), 4–6. https://doi.org/10.1038/nclimate1783\nL89 Regarding the reanalysis, ECCO2 has a horizontal resolution of 0.25ºx0.25º, which\nis coarse for comparison with an ocean model of such high resolution as CESM-H. I\npropose that you perform a comparison with the GLORYS model from Copernicus,\nwhich has a resolution of 0.083ºx0.083º and may better align with the model.\nThanks for your comment. The primary reason for utilizing ECCO2 for comparison is\nbecause its direct output of the vertical velocity (w) variable. While GLORYS has a\nhigher spatial resolution, it lacks the w variable output. Nevertheless, we calculated w\nusing the horizontal velocity derived from GLORYS via the continuity equation 𝑤(𝑧 =\n0 𝜕𝑢 𝜕𝑣\n−ℎ) = −∫ ( + )𝑑𝑧 (Tarry et al., 2022). The results demonstrate agreement\n−ℎ 𝜕𝑥 𝜕𝑦\nwith CESM-H both in pattern and magnitude (Fig. R3). The computation of the\ncontinuity equation, however, is sensitive to the adopted deviational method, which\nmay introduce certain uncertainties. Therefore, despite the higher spatial resolution of\nGLORYS, we choose ECCO2 as the validation products due to its direct output of w.\nFigure R3. Same as Fig. 1c-f but for results derived from GLORYS.\nTarry, D. R., Ruiz, S., Johnston, T. M. S., Poulain, P.-M., Özgökmen, T., Centurioni, L.\nR., et al. (2022). Drifter Observations Reveal Intense Vertical Velocity in a Surface\nOcean Front. Geophysical Research Letters, 49(18), e2022GL098969.\nhttps://doi.org/10.1029/2022GL098969\nL87-91 If I am not mistaken, the model overestimates the UI compared to the reanalysis\ndata (supplementary Fig. 1 and Fig. 1). Please add a brief sentence discussing this fact.\nAlthough the differences are within an acceptable range, it is important to note in the\nmanuscript to encourage the modelling community to strive for improvements in the\nfuture.\nThanks for your advices. We fully agree with you that the simulated UI is overestimated\nby the CESM-H compared with the ECCO2 reanalysis. As noted in previous studies,\nthe coastal upwelling zone simulated in oceanic models may be set by the parameterized\nphysical or numerical horizontal mixing via their effects on the horizontal Ekman layer\nwidth. Further increasing the oceanic resolution would probably make the coastal\nupwelling stronger and narrower (Chang et al., 2023). Thus, the upwelling intensity in\nECCO2 may be underestimated due to the coarse resolution of its oceanic model.\nIndeed, in GLORYS, which has a higher resolution, the simulated upwelling intensity\nis more comparable to that in the CESM-H (Fig. R3), leading further supports to our\nhypothesis. However, even CESM-H is not fine enough to well resolve the coastal\nupwelling. We have cautioned readers about this issue. Please see lines 104-106.\nChang, P., Xu, G., Kurian, J., Small, R. J., Danabasoglu, G., Yeager, S., et al. (2023).\nUncertain future of sustainable fisheries environment in eastern boundary\nupwelling zones under climate change. Communications Earth & Environment,\n4(1), 1–9. https://doi.org/10.1038/s43247-023-00681-0\nL94 I find the methodology proposed by the authors very interesting. However, I have\na question about how you calculate the onset and termination of upwelling. Typically,\nthe change from positive to negative in upwelling sign can be influenced by local wind\npatterns or even mesoscale processes from one day to the next. Have the authors applied\nany filters or analyses to calculate the average sign change over several days? Or have\nyou simply calculated the sign change from one day to the next, and this does not vary\nfrom onset to termination?\nAs you said, the upwelling index (UI) sign could be influenced by local wind patterns\nor even mesoscale processes from one day to the next. Thus, it is difficult to detect the\nprecise transition day from the downwelling to the upwelling. Here, we have applied a\n1-year 4th-order Savitzky-Golay filter (Savitzky & Golay, 1964) to remove the inter-\nseason signals of UI before the detection of onset and termination dates. This would\nhelp us to smooth the seasonal cycle and avoid the high-frequency fluctuation\ninfluences. In the revised manuscript, we have added this information to make the\nmethodology clearer. Please see Line 462-465.\nSavitzky, Abraham., & Golay, M. J. E. (1964). Smoothing and Differentiation of Data\nby Simplified Least Squares Procedures. Analytical Chemistry, 36(8), 1627–1639.\nhttps://doi.org/10.1021/ac60214a047\nL95 I have a question about the choice of poleward regions. In this analysis, the\nBenguela upwelling does not appear as the authors have not subdivided this system.\nWhy? The BCS has the same latitudinal extent as the California upwelling, and the\nseasonal latitude-dependent difference is distinguishable in Fig. 1f. For example, it\ncould be divided at 27ºS. In this regard, please also include the BCS in the rest of the\nanalyses, such as in Fig. 2.\nThanks for your comment. In the revised manuscript, we have subdivided the Benguela\nsystem into three regions (poleward (-P), central (-C), and equatorward (-E)) based on\ndifferences in upwelling seasonality, intensity as well as its potential dynamics. In\nparticular, the BCS-E (15°S-22°S) is located at the poleward part of Benguela system,\nexhibiting a significant seasonal variability, with a greater magnitude during the winter\nand a lesser magnitude during the summer, respectively. The BCS-C region (22°S-29°S)\nexhibit a strong permanent upwelling but weak seasonal variability. As for the BCS-P\nregion (29°S-34°S), situated at the southernmost tip of the African continent, also\nexhibits significant seasonal variability. However, the upwelling seasonal variability\nexhibits a greater magnitude during summer and a lesser magnitude during winter,\nwhich is in contrast to that in BCS-E.\nFollowing the new subdivision method, we have included BCS in the rest of the\nanalyses. However, different from other systems, the BCS region exhibits a permanent\nupwelling in every subdomain. Furthermore, the permanent upwelling in BCS is\nprojected to continue through the 21st century, indicating non-anticipated changes in\nthe timing and duration of the upwelling season. Therefore, we have subtracted the BCS\nregion from the timing analysis as shown in Fig 2.\nL95. Add information to the graph to clarify that the y-axis of Fig. 2a-c represents dates.\nThanks for your suggestion, we have added the information in the revised Fig. 2 and its\ncaption.\nL117. The study is conducted with the CESM-H model forced by the RCP8.5 scenario\nfrom CMIP5, but the comparison with global models is from CMIP6 (it is not\nmentioned in the text whether it is under the SSP5-8.5 scenario). Are there differences\nbetween the different scenarios? Add a sentence to the text or methodology discussing\nthis.\nThanks for your comment. As you said, the CESM-H model is forced by the RCP8.5\nscenario, while the CMIP6 is forced by the SSP5-8.5 scenario. Compared with the\nRCP8.5 scenario, the SSP5-8.5 scenario is also combined with the socioeconomic\nreasons. However, considerable consistency was observed between these two scenarios\n(Gulev et al., 2021), enhancing confidence in comparing the results from CESM-H with\nCMIP6. We have added this discussion in the revised manuscript, please see in line\n402-407.\nL133 Change Pacific to Atlantic\nRevised.\nL140-146 Please add the Benguela subdivision to compare whether it also increases, as\nin the Chile region.\nThanks for your comments. Following your advice, we have subdivided the Benguela\nsystem into three subdomains (See our reply to your comment L95). As shown in Fig.\nR4, all the subdomains in the Benguela are projected to experience an attenuation in\nupwelling intensity. However, the magnitude of the decreasing trends is highly\ndependent on latitudes. In particular, the equatorward region in BCS (BCS-E) is\nprojected to experience a significant attenuation in upwelling intensity at a rate -12.4%\nper century. As for the central region in BCS (BCS-C), the upwelling intensity also\nexhibit a decreasing trend (-7.1% per century), but with a weaker magnitude than that\nin the BCS-E. In contrast, the BCS-P, located at the poleward portion of the South\nAtlantic, exhibits a very faint negative trend (-1.0% per century) in upwelling intensity.\nThis decreasing trend with increasing latitudes is consisted with that in the poleward\nand central regions of HCS located at similar latitudes.\nDifferent from other upwelling systems, the BCS region is located at a lower latitude\ndue to the geographical limitation of African continent. Thus, the increasing intensity\ntrend in the Chile region (renamed as HCS-P in the revised manuscript) could not be\ndetected in the BCS region, due to its location in a higher latitude. However, based on\nthe observed pattern of decreasing upwelling intensity in the BCS with latitudes, it can\nbe inferred that if the Africa continent extents to higher latitudes, there is a potential for\na future increasing trend in upwelling intensity at higher latitudes.\nThis discussion has been added in the revised manuscript, please see in line 197-198.\nFig. R4 Projected changes in upwelling intensity in EBUSs.\nL169-173 I do not believe this study provides results to claim that it can assess local\nnutrient effects. By averaging over hundreds of kilometers both latitudinally and\nlongitudinally, the ability to evaluate local effects is compromised (see Lovecchio et al.\n2017, 2018; Santana-Falcón et al. 2020). Please avoid mentioning that local changes in\nnutrients are being assessed.\nLovecchio E, Gruber N, Münnich M, Lachkar Z (2017) On the longrange offshore\ntransport of organic carbon from the Canary Upwelling System to the open North\nAtlantic. Biogeosciences 14(13):3337–3369. https:// doi. org/ 10. 3929/ ethz-b- 00019\noffshore transport of organic carbon from the Canary Upwelling System to the open\nNorth Atlantic. Biogeosciences 15:5061–5091. https:// doi. org/ 10. 5194/ bg- 15- 5061-\nupwelling filaments in the Canary Current System. Prog Oceanogr\n184(April):102322. https://doi.org/10.1016/j.pocean.2020.102322\nSorry for our misleading statement. We have revised this statement in the revised\nmanuscript to avoid misleading readers.\nL173-174 I understand that the CESM-H model lacks biogeochemistry, so the authors\nuse the global CMIP6 models to test the hypothesis. Although the results obtained\npreviously by CMIP6 models could be compared with those of CESM-H, no evaluation\nof the models has been provided to see if they resolve latitudinal upwelling as in the\ncase of CESM-H. In fact, there are different studies (Varela et al. 2022; Sylla et al. 2022)\nshowing how CMIP6 models differ from observations. Therefore, I believe this section\nprovides useful information that complements previous results and helps take a first\nstep regarding the biogeochemical aspect of upwelling. However, a paragraph should\nbe added discussing the deficiencies of CMIP6 models in reproducing EBUSs,\nprimarily associated with the lack of resolution.\nThanks for your comments and useful references. As evidenced by previous studies, the\naccuracy of upwelling simulation is contingent upon on realistic representation of the\ncross-shore wind stress structure, seafloor topography, and oceanic mesoscale eddies,\nwhich all connect to the high-resolution in oceanic or atmospheric model configuration.\nThus, using the low-resolution CMIP6 models would entail large uncertainties in\nprojecting the future changes in upwelling intensity and associated biogeochemical\nchanges. In the revised manuscript, we have cautioned readers about this issue. Please\nsee Lines 303-306.\nL185-186 The mentioned section in the methods does not exist.\nRevised.\nL216-219 Do you compare the Mixed Layer Depth (MLD) and nutrients based on\naverages from the areas shown in Fig 1? If so, I believe changes in MLD may be\ninadequately represented as it can be highly variable from one grid cell to another.\nAdditionally, the authors do not assess stratification at any point. I suggest they\ncalculate the Brunt-Väisäla frequency as an indicator of stratification changes (see\nVázquez et al. 2023)\nVázquez, R., Parras-Berrocal, I. M., Koseki, S., Cabos, W., Sein, D. V., & Izquierdo, A.\n(2023). Seasonality of coastal upwelling trends in the Mauritania-Senegalese region\nunder RCP8. 5 climate change scenario. Science of the Total Environment, 898, 166391.\nThanks for your comment and useful reference. Following the method in Vázquez et al.\n(2023), we have calculated the upper 150 m Brunt-Väisäla frequency as an indicator of\nstratification changes instead of the MLD in the revised manuscript. As shown in Fig\nR5, the stratification in major EBUSs exhibits pronounced increasing trends. However,\nlimited consensus was also exhibited regarding the relationship between Brunt-Väisäla\nfrequency changes and primary production changes. Therefore, we can still conclude\nthat the local upwelling changes may overwhelm stratification changes, governing the\necosystem’s productivity. Moreover, we also performed sensitivity analysis to evaluate\nthe influence of choosing different average depths in calculating the stratification (not\nshown). They all show less relationship between stratification changes and primary\nproduction changes, leading supports to the validity of the conclusions in our\nmanuscript. The above analysis has added in the revised manuscript, see lines 292-296.\nFig. R5 Same as Fig. 4c, but for the relationship between Brunt-Väisäla frequency\nchanges and primary production changes.\nMajor comments:\nThe first question I have regarding the results in this work is related to the use of model\nvertical velocity outputs, as well as the use of an arbitrary depth (50 m) as the most\nrepresentative of upwelling. I would like the authors to justify why they chose to use\n50 m. In studies on upwelling systems, it has been agreed upon that there is significant\nvariability both latitudinally and seasonally in the upwelling source depth. Using 50 m\nas a criterion oversimplifies the complexity of EBUSs. In this regard, using the MLD\nas the depth at each latitude and for each season would be more appropriate. However,\nestimating the appropriate depth remains highly complex. In fact, Jacox et al. (2018)\nproposes another method to assess upwelling volume that is superior to model vertical\nvelocity outputs (Marchesiello and Estrade, 2010; Rossi et al. 2013; Oerder et al. 2015).\nIn this context, I include the paragraph that describes this issue:\n“However, we do not use modeled vertical velocities to construct the index for several\nreasons: First, the calculation of upwelling using modeled vertical velocity is not\nstraightforward. One must first define a representative depth at which vertical velocity\nis extracted (e.g., the MLD), which will be different for each grid cell, making it\ndifficult to close the transport budget. For example, when the MLD differs between\nadjacent grid cells, water can enter/exit the surface mixed layer horizontally and that\ncomponent of the transport will be missed from the upwelling/downwelling estimate.”\nTherefore, I would like to suggest to the authors that they reconsider the use of the MLD\nas the depth criterion and discuss the use of vertical velocity versus the wind field.\nMarchesiello, P., and P. Estrade (2010), Upwelling limitation by onshore geostrophic\nflow, J. Mar. Res., 68, 37–62\nRossi, V., Feng, M., Pattiaratchi, C., Roughan, M., & Waite, A. M. (2013). On the\nfactors influencing the development of sporadic upwelling in the Leeuwin Current\nsystem. Journal of Geophysical Research: Oceans, 118(7), 3608-3621.\nOerder, V., Colas, F., Echevin, V., Codron, F., Tam, J., Belmadani, A., 2015. Peru-Chile\nupwelling dynamics under climate change. J. Geophys. Res. Oceans 120 (2), 1152–\n1172. https://doi.org/10.1002/2014JC0102, 99.\nJacox, M.G., Edwards, C.A., Hazen, E.L., Bograd, S.J., 2018. Coastal upwelling\nrevisited: Ekman, Bakun, and improved upwelling indices for the U.S. West Coast. J.\nGeophys. Res. Oceans 123, 7332–7350. https://doi.org/10.1029/2018JC014187.\nThanks for your comments and useful references. We fully agree with you that choosing\nmixed layer depth (MLD) as the representative depth of vertical velocity is more\nappropriate for studying the upward transport of water due to the upwelling. However,\nwe choose 50 m as the representation of upwelling in this study for the following two\nreasons. First, due to storage limitation, daily averaged oceanic vertical velocity is only\noutput at selected depths (50, 105, 528 and 1146 m) in CESM-H. Considering that the\nMLD at the EBUSs approximates 50 m, we have therefore selected this value as the\nrepresentation depth for upwelling in our study. Second, as you pointed out, the MLD\nvaries with latitude and seasons in the EBUSs, posing challenges in closing the\ntransport budget due to the discontinuous distribution of MLD. Therefore, we have\nopted to analyze the future changes in upwelling characteristics, encompassing its\ntiming and intensity, at a fixed depth of 50 m, which closely approximates the MLD\nwithin the EBUSs.\nHowever, we have performed a sensitivity analysis to evaluate the influence of\nchoosing different depths in projecting future changes in upwelling characteristics. It is\nworth noting that the CESM-H only provides daily vertical velocity data at a depth of\n50 m, therefore, we utilized monthly data for conducting the sensitivity analysis.\nAccording to observation and our CESM-H simulation, the MLD in EBUSs ranges\nfrom 30 m to 60 m. Figure R6 displays the climatological upwelling index as well as\nits changes in the future at 30m, 50m, and 60m, respectively. Although some\nquantitative variations are observed among the upwelling changes at different depths,\nthey consistently demonstrate an increase in UI during the month of upwelling,\nindicating a progressive advancement of the onset date. Similar UI changes among\ndifferent depths is also detected at the termination month of upwelling across all the\nEBUSs. Moreover, the earlier onset and prolonged duration of upwelling season are\nconsistently observed in CMIP6 regardless of the chosen depth (Fig. R7). The\nupwelling intensity changes at different depths are also compared, yielding highly\nconsistent results that robustly support the validity of our manuscript's conclusions (Fig.\nR8).\nFig. R6. Climatology UI (black line) and future UI trend near transition months (orange\nbar) at 30 m (a-d), 50 m (e-h), and 60 m (i-l).\nFig. R7, Same as Fig. 2e-h, but for 30 m (a-d), 50 m (e-h), and 60 m (i-l).\nFig. R8 Same as Fig. 3e, but for 30 m (a), 50 m (b), and 60 m (c).\nWe opt for utilizing modeled vertical velocities instead of the wind-derived upwelling\nindex due to the following reasons. First, the wind-derived upwelling index (UI ) is\ne\ndepended on the assumption that coastal upwelling is majorly induced by the wind\nstress. This assumption could hold when the time scale is shorter, as proved in\nobservations (Chereskin, 1995). However, over longer time scales, oceanic and\natmospheric processes beyond wind stress may exert a significant influence on the\nupwelling variation (Ding et al., 2023; Jing et al., 2023; Rykaczewski et al., 2015). In\nparticular, recent studies have demonstrated that oceanic processes may surpass the\ninfluence of wind stress, playing a pivotal role in modulating long-term changes in\nupwelling under greenhouse warming (Jing et al., 2023). Second, a key motivation of\nthis study lies in the assessment of shifting rates of upwelling onset and termination\nunder greenhouse warming. However, our findings demonstrate significant disparities\nin estimating the upwelling onset and termination between the modeled vertical velocity\nand wind-derived upwelling index (Fig R9). Thus, employing the wind-derived\nupwelling index for assessing future changes in upwelling timing and intensity may\nintroduce considerable uncertainties.\nFig. R9 Climatology UI derived from vertical velocity (black) and wind (red).\nIn the revised manuscript, we have added the sensitive analysis in selecting\nrepresentative depths for assessing changes in upwelling timing and intensity.\nAdditionally, we have emphasized the limitations of using a fixed depth instead of MLD\nwhen estimating upwelling, see line 451-458.\nChereskin, T. K. (1995). Direct evidence for an Ekman balance in the California Current.\nJournal of Geophysical Research: Oceans, 100(C9), 18261–18269.\nhttps://doi.org/10.1029/95JC02182\nDing, H., Alexander, M. A., & Ting, M. (2023). Revisiting the Relationship between\nthe North Pacific High and Upwelling Winds along the West Coast of North\nAmerica in the Present and Future Climate. Journal of Climate, 36(23), 8211–\n8224. https://doi.org/10.1175/JCLI-D-23-0238.1\nJing, Z., Wang, S., Wu, L., Wang, H., Zhou, S., Sun, B., et al. (2023). Geostrophic flows\ncontrol future changes of oceanic eastern boundary upwelling. Nature Climate\nChange, 1–7. https://doi.org/10.1038/s41558-022-01588-y\nRykaczewski, R. R., Dunne, J. P., Sydeman, W. J., García-Reyes, M., Black, B. A., &\nBograd, S. J. (2015). Poleward displacement of coastal upwelling-favorable winds\nin the ocean’s eastern boundary currents through the 21st century. Geophysical\nResearch Letters, 42(15), 6424–6431. https://doi.org/10.1002/2015GL064694\nThe second concern is related to the domain chosen for the zonal average. In this regard,\nthe authors average over 200 km offshore and mention that this is most appropriate for\ncapturing both Ekman transport-associated upwelling and upwelling associated with\nEkman transport divergence. Averaging over such a wide portion of the coast may\nobscure the upwelling transport. Perhaps with low-resolution global models, obtaining\na 200 km average is necessary. However, with a model of such high resolution as\nproposed in this study, the authors should consider narrowing this region to a maximum\nof 100 km. This would allow for a more appropriate assessment of upwelling, where\nmesoscale processes play a significant role. Therefore, I propose that the authors\nconduct the analysis by limiting the average to the first 100 km of the coast. In fact, in\nthe methodology, the authors mention the work of Du et al. (2023), who also conducted\na study with an average over the first 50 km offshore. Moreover, in Fig. 1a of the\nmanuscript, the vertical velocity associated with upwelling is closely confined to the\ncoast.\nThanks for your comments. As you pointed out, the simulated upwelling in the four\nEBUSs by CESM-H consists of a rapid coastal component driven by the alongshore\nwind stress and confined to a narrow band (~50 km) next to the coast and a slower\ncomponent driven by the wind stress curl and extending further off shore (Fig 1a),\nwhich is consistent with the existing theoretical arguments and high-resolution regional\nocean simulations. The simulated upwelling by the CMIP6 models, however, exhibits\na weaker and more diffusive nature due to its coarser resolution, resulting in the absence\nof the prominent coastal upwelling zone. Thus, integrating vertical velocity over a\nsufficiently wide region should be less sensitive to model resolution and more suitable\nfor the CMIP6 analysis. For this reason, we integrate the vertical velocity within ~200\nkm from the coast in the individual EBUSs and evaluate their long-term changes under\ngreenhouse warming. Such integration covers the coastal upwelling and a large fraction\nof offshore wind stress curl-driven upwelling (Jacox et al., 2018), both of which are\nsuggested to influence the ecosystem.\nFollowing your advice, in the revised manuscript, we recompute the UI for the\napproximately 50-km-wide coastal zone in order to reveal the response of coastal\nupwelling to greenhouse warming. The secular changes in the timing and intensity of\nupwelling in the individual EBUSs resemble those calculated from the approximately\n200-km-wide upwelling band from the coast (Fig. R10 and R11). The coastal upwelling\nwithin the 50-km-wide coastal zone from the coast also exhibits a seasonal\nadvancement and prolonged duration, mirroring those found in the 200-km-wide\ncoastal zone from the coast. Regarding the upwelling intensity, its secular changes also\nexhibit complicated variations across the EBUSs. However, the secular changes in\nindividual EBUSs bear resemblance to those calculated at a distance of 200 km from\nthe coastline, albeit with a greater magnitude. Such resemblance suggests that the\nresponse of upwelling timing and intensity under greenhouse warming is a robust\nfeature that is not sensitive to the width of the upwelling band selected for analysis. We\nhave added the above discussion in the revised manuscript, please see in line 128-134\nand line 201-206.\nFig. R10. Same as Fig. 2a-d, but for results averaged within 50 km offshore.\nFig. R11. Same as Fig. 3a-d, but for results averaged within 50 km offshore.\nJacox, M. G., Edwards, C. A., Hazen, E. L., & Bograd, S. J. (2018). Coastal Upwelling\nRevisited: Ekman, Bakun, and Improved Upwelling Indices for the U.S. West\nCoast. Journal of Geophysical Research: Oceans, 123(10), 7332–7350.\nhttps://doi.org/10.1029/2018JC014187\nReply to the second reviewer\nWe are very grateful to you for your time in carefully reading our manuscript and\nproviding helpful comments that make our manuscript better. We have carefully\nconsidered each of your comments (in blue) and revised the manuscript accordingly.\nPlease find our response (in black) to your comments below.\nReviewer #2 (Remarks to the Author):\nThis paper is a model-based study (CESM-H, CMIP6, ECCO2) aimed at understanding\ncoastal upwelling variability in the context of the climate change. While the authors\nhave undertaken significant efforts to comprehend the seasonal progression and\nupwelling intensity in the EBUSs, there are concerns regarding certain processes and\nmodel biases than prevent me for approving its publication.\nMajor concerts:\n1) Recently, it has been pointed out that the geostrophic flow may play a significant\nrole in controlling the upwelling trends under the climate change (Jing et al., 2023). The\nauthors should include an explanation regarding whether and why it is not being\nconsidered in this study. Could the results be the same if the geostrophic flow is\nconsidered?\nThanks for your comments. In this study, we have utilized the modelled vertical velocity\n(w) outputs to represent the coastal upwelling processes, which allows us to account for\nthe contribution of geostrophic flow, as indicated in previous studies (Ding et al., 2021;\nDu et al., 2023; Jing et al., 2023). In accordance with your guidance, we have conducted\nfurther analyses to disentangle the impacts of wind and geostrophic flow on future\nupwelling changes based on CESM-H, thereby elucidating the underlying mechanisms\ndriving future upwelling changes under greenhouse warming. It is noteworthy that the\nanalysis presented herein pertains solely to the monthly data due to the unavailability\nof daily averaged three-dimensional temperature and salinity. However, we were able\nto assess the upwelling timing changes by examining the variations in upwelling during\nthe transitional month: a positive linear trend of upwelling indicates an earlier onset or\na delay in termination, while a negative trend suggests a later onset or an advancement\nin termination. As show in Fig R1, there is a uniform increasing trend in the upwelling\nduring upwelling onset month across all EBUSs, indicating an earlier-shifting trend of\nupwelling occurrence. However, the primary contributor varies across different ocean\nbasins. In the Pacific basin, the positive trend of upwelling during the onset month is\nprimarily attributed to changes in wind-driven Ekman transport (UI ), illustrating the\ne\ndominate role of wind variations in driving the advancement of upwelling season. In\ncontrast, in the Atlantic basin, the geostrophic transport (UI ) emerges as the primary\ng\ncontributing factor to the increasing trend in the positive trend of upwelling during the\nonset month. Similar situation is also detected for the upwelling changes during the\nterminal month of upwelling season across all the EBUSs.\nFig. R1 The climatology UI (black) and UI trend during 1920-2100 (purple) near\ntransition months. The wind-induced UI (UI ) trend and geostrophic-induced UI (UI )\ne g\ntrend were shown in red and blue. The orange tringle demonstrates the onset of\nupwelling while the inverted triangle represents the termination date.\nIn addition, the same decomposition is also conducted for the future changes in the\nupwelling intensity (Fig. R2). Consistent with a recent study (Jing et al. 2023), in most\nEBUSs, the secular changes in upwelling intensity is primarily attributed to the\ngeostrophic flow changes, highlighting the controlling role of geostrophic flows in\nupwelling intensity trends in EBUSs under greenhouse warming. However, as latitudes\nincrease, the dominant influence of geostrophic flow gradually diminishes, yielding to\nthe prominence of wind-induced processes. Specifically, in the poleward region of\nCanary Current system and Humboldt Current system, the wind-derived Ekman\ntransport exceeds the geostrophic flow, determining the future changes in upwelling\nintensity.\nFig. R2. Trend of upwelling intensity calculated by UI (purple), UI (blue), and UI (red)\ng e\nderived from CESM-H. The error bar represents the standard error of linear regression.\nWe have the above analysis and discussion in the revised manuscript. See lines 156-\n175 and 233-243.\nDing, H., Alexander, M. A., & Jacox, M. G. (2021). Role of Geostrophic Currents in\nFuture Changes of Coastal Upwelling in the California Current System.\nGeophysical Research Letters, 48(3), e2020GL090768.\nhttps://doi.org/10.1029/2020GL090768\nDu, T., Wang, S., Jing, Z., Wang, H., & Wu, L. (2023). Greenhouse Warming Weakens\nthe Seasonal Cycle of the Eastern Boundary Upwelling. Geophysical Research\nLetters, 50(11), e2023GL103857. https://doi.org/10.1029/2023GL103857\nJing, Z., Wang, S., Wu, L., Wang, H., Zhou, S., Sun, B., et al. (2023). Geostrophic flows\ncontrol future changes of oceanic eastern boundary upwelling. Nature Climate\nChange, 1–7. https://doi.org/10.1038/s41558-022-01588-y\n2) For the computation of the upwelling timing and intensity, the vertical velocity at 50\nm depth has been considered. Could the results be different using another depth? (Note\nthat Xiu et al. (2018) shows that, at least for CCS, there are three major drivers of the\nvertical velocity changes which are depth-depending: alongshore winds, wind-curl and\nmesoscale activity)\nFollowing your comments, we have performed a sensitivity analysis to evaluate the\ninfluence of choosing different depths in projecting future changes in upwelling\ncharacteristics. It is worth noting that the CESM-H only provides daily vertical velocity\ndata at a depth of 50 m; therefore, we utilized monthly data for conducting the\nsensitivity analysis. According to observation and our CESM-H simulation, the mixed\nlayer depth in EBUSs ranges from 30 m to 60 m. Figure R3 displays the climatological\nupwelling index as well as its changes in the future at 30m, 50m, and 60m, respectively.\nAlthough some quantitative variations are observed among the upwelling changes at\ndifferent depths, they consistently demonstrate an increase in UI during the month of\nupwelling, indicating a progressive advancement of the onset date. Similar UI changes\namong different depths is also detected at the termination month of upwelling across\nall the EBUSs. Moreover, the earlier onset and prolonged duration of upwelling season\nare consistently observed in CMIP6 regardless of the chosen depth (Fig. R4). The\nupwelling intensity changes at different depths are also compared, yielding highly\nconsistent results that robustly support the validity of our manuscript's conclusions (Fig.\nR5). These has been added in line 451-458.\nFig. R3. Climatology UI (black line) and future UI trend near transition months (orange\nbar) at 30 m (a-d), 50 m (e-h), and 60 m (i-l).\nFig. R4, Same as Fig. 2e-h, but for 30 m (a-d), 50 m (e-h), and 60 m (i-l).\nFig. R5 Same as Fig. 3e, but for 30 m (a), 50 m (b), and 60 m (c).\n3) How do the VNT and NPP appear under historical period compared to observational\ndata (e.g., 2003-2022)? Most of the models exhibit a significant bias in the EBUSs,\nparticularly with the biogeochemical variables. Can we have confidence in the results\nconsidering these biases?\nThanks for your comments. The sparsity of oceanic observations on ocean vertical\nvelocity measurements leads to difficult in measure VNT directly. Therefore, we only\nevaluate the simulated performance of NPP by ESMs compared to the observations (Fig.\nR6). We agree with you that the ESMs in CMIP6 is not fine enough to well simulate\nthe NPP. In particular, the simulated NPP in the EBUSs is weaker and more diffusive\ncompared with the observation. Significantly, the simulation bias in NPP bears a\nstriking resemblance to that observed in low-resolution modes for coastal upwelling\n(Fig. 1 in Jing et al. 2023), suggesting that the NPP bias in ESMs may primarily stem\nfrom inadequate representation of coastal upwelling in these models. Further increasing\nthe oceanic resolution would probably make the coastal upwelling stronger and\nnarrower, thereby mitigating the simulation bias in NPP within the EBUSs. In addition,\nthe low-resolution ESMs used in this study cannot resolve the oceanic small-scale\nprocesses (e.g., eddy activities, coastal trapped waves). Although resolving these small-\nscale processes remains probably unaccomplishable for the long-term global climate\nsimulations up to now, they are suggested to contribute notably to the vertical nutrition\ntransport in the mixed layer and play an important role in modulating the local NPP.\nHence, the simulation bias of coastal upwelling and the deficiency in accounting for\noceanic small-scale processes would both contribute to uncertainties in projecting\nfuture changes in NPP within EBUSs. We have cautioned readers about this issue in the\nrevised manuscript. Please see Lines 303-306\nFig. R6 a, Observed vertically integrated primary production during 2003-2022, b,\nvertically integrated NPP derived from CMIP6 ESMs.\nMinor points:\nLine 111-114: Is this valid for Southern HCS and Northern HCS region? (as far as it is\nknown, these two regions have different upwelling regimes, which one being permanent\nand the other seasonal)\nFollowing your comment, we have carefully reconsidered the selection of HCS region.\nIn the revised manuscript, the HCS region is defined from 5ºS to 44ºS, and further\nsubdivided into three regions (poleward region (-P), central region (-C), and\nequatorward region (-E)) based on differences in upwelling seasonality, intensity as\nwell as its potential dynamics. The HCS-E (5°S-15°S) is located at the equatorward part\nof the HCS, which exhibits a permanent upwelling with significant seasonal variability.\nSpecifically, it experiences a higher magnitude during winter and a lower magnitude\nduring summer. The HCS-C (15°S-36°S) situated at the central region of the HCS, also\nexhibits a strong permanent upwelling but weak seasonal variability. As for the HCS-P\nregion (36°S-44°S), located at the poleward extent of the HCS, exhibits a seasonal\nupwelling phenomenon, with upwelling events predominantly occurring during\nsummer and downwelling events observed in winter.\nWe have conducted a reanalysis on the future changes in the timing of the upwelling\nseason in the subdivided regions. There is an obvious seasonal advancement and\nextension of coastal upwelling in the HCS-P. As for the HCS-E and HCS-C regions, the\npersistent annual upwelling is projected to remain throughout the 21st century,\nindicating no anticipated changes in the timing and duration of the upwelling season\n(Fig R7).\nFig. R7 The upwelling timing changes in HCS-P, HCS-C, and HCS-E\nLine 235 Fig.1: why are figures 1c and 1d composed of two subfigures? (Should they\nbe like 1f?)\nThis has been revised in the manuscript since we change the definition of different\nsubregions.\nLine 515 Supp. Fig.2b and 2c: the curves in these figures represent the UI only for one\nof the initially selected subdivisions for each upwelling system (as shown in figure 1)?\nSorry for the misleading. The UI curves shown in Supplementary Fig. 2 represent for\nonly one subregion for each EBUSs. We have rewritten the figure caption to make it\nless ambiguous.\nLine 547 Supp. Table 1: Include units on the table\nAdded.\nReply to the third reviewer\nWe sincerely appreciate your thorough review of our manuscript and the insightful\ncomments you have provided. Your feedback has greatly contributed to improving the\nquality of our work. We have thoughtfully addressed each of your suggestions (in blue)\nand incorporated revisions accordingly. Below, you will find our responses (in black)\nto your comments.\nReviewer #3 (Remarks to the Author):\nTitle: Review of ‘Future changes in coastal upwelling and biological production in\neastern boundary upwelling systems’\nThe paper investigates future changes in the duration and intensity of upwelling in the\neastern boundary regions and its potential consequences on primary production.\nUtilizing simulations from CMIP6 models and high-resolution CESM, the authors\nshowed that the duration of upwelling is projected to be longer under high-emission\nscenario. But changes in the upwelling intensity show region-specific features. The\nauthors further explore the implications of these intensity changes on net primary\nproduction. While the results are interesting and seem to be scientifically sound, this\npaper still needs substantial work before it can be published. I summarize my concerns\nand comments below.\nMajor concerns:\n1) Refinement of introduction. This paper introduces a novel investigation into changes\nin upwelling duration, a distinctive aspect that warrants attention. However, it is crucial\nto incorporate relevant findings from existing literature to enrich the contextual\nbackground. For instance, Wang et al. (2015), utilizing CMIP5 models, observed an\nincreased duration of upwelling in the high latitudes of major EBUSs, with differing\ntrends in low latitudes. Summarizing these findings explicitly in the introduction is\nessential to provide readers with a comprehensive understanding of the advancements\nin this field. In addition, the introduction of future changes in net primary production\nonly discussed ocean stratification and upwelling intensity. However, nutrients in the\nsource water can also impact future changes in net primary production (Rykaczewski\nand Dunne, 2010). Please make the introduction more comprehensive.\nThanks for your comments and useful references, we have refine the introduction to a\nmore comprehensive version in line 42-70.\n“Consequently, Bakun hypothesized that intensified land-sea thermal contrast under\nfuture warming condition would produce stronger pressure gradients, thereby\nenhancing coastal upwelling in the EBUSs by fortifying alongshore winds (Bakun,\n1990). However, recent studies have underscored the impact of poleward migration of\nthe high-pressure systems on future changes in coastal upwelling, projecting an\nintensified upwelling in higher latitudes and weakened upwelling in lower latitudes\n(Rykaczewski et al., 2015). Besides, under greenhouse warming, an earlier and\nprolonged upwelling season was also projected in the high latitudes of EBUSs, while a\ndifferent trend was identified in the low latitudes (Wang et al., 2015). These studies\nequate changes in upwelling with changes in the upwelling-favorable winds. However,\na meta-analysis revealed substantial diversity in the responses of regional coastal\nupwelling to anthropogenic forcing, even under universally increasing wind trends in\nmost EBUSs (Abrahams et al., 2021; Brady et al., 2017; Jing et al., 2023; Oyarzún &\nBrierley, 2019; Sousa et al., 2020; Sydeman et al., 2014; Tim et al., 2015). In particular,\nfuture changes in stratification (Oyarzún & Brierley, 2019; Sousa et al., 2020; Vázquez\net al., 2023) and geostrophic transport (Ding et al., 2021; Jing et al., 2023) could surpass\nthe impact of winds, playing an important role in modulating the coastal upwelling.\nTherefore, considerable uncertainty persists regarding the evolution of coastal\nupwelling with climate change in terms of the intensity and timing (García-Reyes et al.,\n2015; Wang et al., 2015).”\n“It is generally thought that greenhouse warming should reduce the productivity in the\nEBUSs by inhibiting the upward nutrient supply due to intensified stratification, which\nis supported by a majority of studies (Fu et al., 2016; McGowan et al., 2003; Oyarzún\n& Brierley, 2019; Roemmich & McGowan, 1995; Sousa et al., 2020). However, a lack\nof consensus on future productivity projections derived from climate models challenges\nthe notion that stratification is the primary determinant in predicting forthcoming\nproductivity changes in the EBUSs (Bograd et al., 2023). Evidence exists that the\ncoastal upwelling contributes nonnegligible to the long-term changes of productivity in\nsome EBUSs (Bakun et al., 2015; Xiu et al., 2018). Besides, potential alterations in the\nnitrate content of upwelling source water could also exert a significant influence on\nproductivity (Renault et al., 2016; Rykaczewski & Dunne, 2010; Xiu et al., 2018).\nHowever, how the evolving dynamics of upwelling affect the nutrient supply and\nproductivity in the major EBUSs remains shrouded in uncertainty.”\n2) Refinement of Chile and Canary-Iberian upwelling systems definition. The definition\nof the Chile and Canary-Iberian upwelling systems requires careful consideration based\non the results presented in Figure 1a. The depicted large upward vertical velocity\nextends beyond the designated boxes in the Chile and Canary regions, reaching\nsubstantially equatorward. Notably, the surface wind in the Chile box is perpendicular\nto the coast, which may not accurately represent the typical Southern HCS upwelling\nregions. Additionally, it is observed that stronger coast-parallel winds (as shown in\nFigure 1a) and higher net primary production (as indicated in Figure 1b) occur south of\nthe defined Canary box.\nIn addition, if authors define upwelling regions based on vertical velocity, Figure 1a\ndemonstrates that positive vertical velocity has very narrow zonal extents. This suggests\nthat selecting 200 km away from the coast will include regions with negative vertical\nvelocity. Given that upwelling timing and intensity calculations rely on box averages,\nthe selection of these boxes has the potential to influence the results of future changes\nin upwelling intensity and duration. To ensure the accuracy of the study's findings,\na reconsideration of the chosen boxes for defining the Chile and Canary-Iberian\nupwelling systems is recommended. This adjustment should aim to better capture the\nrepresentative characteristics of these upwelling regions and minimize potential biases\nin the assessment of future changes in upwelling dynamics.\nThanks for your suggestions, we have expanded our domains according to your advice\nand included the definition into the methodology, please see in line 422-443.\n“In terms of seasonality, each EBUS is divided into three regions: the poleward region\n(-P), the central region (-C), and the equatorward region (-E). The poleward regions,\nlocated at high latitudes in each EBUS, experience a pronounced seasonal upwelling,\nwith stronger upwelling predominantly occurring during summer and weaker upwelling\nor downwelling observed in winter. The equatorward regions, situated at low latitudes\nin each EBUS, also exhibit a significant seasonal upwelling. However, the upwelling\nreaches its peak during summer and weakens or even shifts to downwelling during\nwinter, in contrast to the poleward regions. The central regions, positioned between the\npoleward and equatorward regions, exhibit a strong permanent upwelling but display\nlimited seasonal variability across all EBUSs.\nAccording to the aforementioned categorization principle, the CanCS is divided into\nthree regions: the poleward region (CanCS-P), situated along the coast of the Iberian\nPeninsula (36°N-43°N), the central region (CanCS-C), encompassing latitudes ranging\nfrom 19°N-36°N, and the equatorward region (CanCS-E), traditionally referred to as\nthe Mauritania-Senegalese upwelling region, spanning from 12°N-19°N. Similarly, the\nHCS located in the South Pacific is also divided into three regions: HCS-E (5°S-15°S),\nHCS-C (15°S-36°S), and HCS-P (36°S-44°S). Likewise, the BCS situated in the South\nAtlantic is categorized as BCS-E (15°S-22°S), BCS-C (22°S-29°S), and BCS-P (29°S-\n33°S). The exception lies in the CalCS, where no winter-intensified seasonal upwelling\nwas observed in its equatorward areas. Consequently, the CalCS region is divided solely\ninto two regions: the central region (CalCS-C) spanning from 30°N to 40°N and the\npoleward region (CalCS-P) ranging from 40°N-48°N.”\nFor the selection of offshore distance, we initially utilized a 200-km average as\nnecessary for the CMIP6 low-resolution models. The simulated upwelling by the\nCMIP6 models, exhibits a weaker and more diffusive nature due to its coarser resolution,\nresulting in the absence of the prominent coastal upwelling zone. Thus, integrating\nvertical velocity over a sufficiently wide region should be less sensitive to model\nresolution and more suitable for the CMIP6 analysis. For this reason, we integrate the\nvertical velocity within ~200 km from the coast in the individual EBUSs and evaluate\ntheir long-term changes under greenhouse warming. Such integration covers the coastal\nupwelling and a large fraction of offshore wind stress curl-driven upwelling, both of\nwhich are suggested to influence the ecosystem.\nHowever, in light of your feedback and to underscore the resolution advantages of\nCESM-H, a 200-km average may not be optimal. Therefore, in the revised manuscript,\nwe recompute the UI for the approximately 50-km-wide coastal zone in order to reveal\nthe response of coastal upwelling to greenhouse warming. The secular changes in the\ntiming and intensity of upwelling in the individual EBUSs resemble those calculated\nfrom the approximately 200-km-wide upwelling band from the coast (Fig. R1 and R2).\nThe coastal upwelling within the 50-km-wide coastal zone from the coast also exhibits\na seasonal advancement and prolonged duration, mirroring those found in the 200-km-\nwide coastal zone from the coast. Regarding the upwelling intensity, its secular changes\nalso exhibit complicated variations across the EBUSs. However, the secular changes in\nindividual EBUSs bear resemblance to those calculated at a distance of 200 km from\nthe coastline, albeit with a greater magnitude. Such resemblance suggests that the\nresponse of upwelling timing and intensity under greenhouse warming is a robust\nfeature that is not sensitive to the width of the upwelling band selected for analysis. We\nhave added the above discussion in the revised manuscript, please see in line 128-135\nand line 201-206\nFig. R1. Same as Fig. 2a-d, but for results averaged within 50 km offshore.\nFig. R2. Same as Fig. 3a-d, but for results averaged within 50 km offshore.\n3) Lack of mechanisms in upwelling changes. One notable limitation in the current\nstudy lies in the absence of a detailed exploration of the underlying mechanisms driving\nchanges in the duration and intensity of upwelling. This omission weakens the novelty\nof the study, especially when compared to the work of Wang et al. (2015), who utilized\nCMIP5 models and already delved into discussions on these changes, demonstrating\ngeneral consistency with the results presented in this paper. To strengthen the scientific\ncontribution and novelty of this study, it is imperative to incorporate a thorough\nexamination of the mechanisms influencing the observed changes in upwelling duration\nand intensity.\nAppreciate for your suggestion. In accordance with your guidance, we have conducted\nfurther analyses to disentangle the impacts of wind and geostrophic flow on future\nupwelling changes based on CESM-H, thereby elucidating the underlying mechanisms\ndriving future upwelling changes under greenhouse warming. It is noteworthy that the\nanalysis presented herein pertains solely to the monthly data due to the unavailability\nof daily averaged three-dimensional temperature and salinity. However, we were able\nto assess the upwelling timing changes by examining the variations in upwelling during\nthe transitional month: a positive linear trend of upwelling indicates an earlier onset or\na delay in termination, while a negative trend suggests a later onset or an advancement\nin termination. As show in Fig R3, there is a uniform increasing trend in the upwelling\nduring upwelling onset month across all EBUSs, indicating an earlier-shifting trend of\nupwelling occurrence. However, the primary contributor varies across different ocean\nbasins. In the Pacific basin, the positive trend of upwelling during the onset month is\nprimarily attributed to changes in wind-driven Ekman transport (UI ), illustrating the\ne\ndominate role of wind variations in driving the advancement of upwelling season. In\ncontrast, in the Atlantic basin, the geostrophic transport (UI ) emerges as the primary\ng\ncontributing factor to the increasing trend in the positive trend of upwelling during the\nonset month. Similar situation is also detected for the upwelling changes during the\nterminal month of upwelling season across all the EBUSs.\nFig. R3 The climatology UI (black) and UI trend during 1920-2100 (purple) near\ntransition months. The wind-induced UI (UI ) trend and geostrophic-induced UI (UI )\ne g\ntrend were shown in red and blue. The orange tringle demonstrates the onset of\nupwelling while the inverted triangle represents the termination date.\nAdditionally, the investigation delves into the potential dynamics of wind changes that\nfavor upwelling in the Pacific EBUSs (Fig. R4). Following the methodology employed\nby Rykaczewski et al. (2015), we demonstrate a significant amplification of wind\npatterns during spring, which exhibits a strong correlation with variations in land-sea\npressure gradients (Fig. R4 a,d). However, these pressure gradient changes are not\nsolely induced by alterations in land-sea thermal contrast as Bakun hypothesized\n(Bakun, 1990). Specifically, the correlation between springtime UI and land-sea\ne\ntemperature differences did not achieve statistical significance (Fig. R4 b,e). Instead,\nthe land-sea pressure gradients exhibit a strong correlation with the atmospheric\npressure in the interior ocean, with their inter-model correlations reaching 0.69 and 0.77\nin the CalCS-P and HCS-P regions, respectively (Fig. R4 c,f). This suggests that the\nintensified along-shore wind during the spring season in the Pacific EBUSs might\nattributed to the large-scale atmospheric circulation changes, which potentially\noriginates from north-south shifts of the ITCZ under a warming climate (Song et al.,\n2018).\nFig. R4 Inter-model relationship between springtime UI and land-sea pressure\ne\ngradients (a,d), land-sea temperature differences (b,e), and atmospheric pressure in the\ninterior ocean (c,f) in CalCS-P and HCS-P.\nIn addition, the same decomposition is also conducted for the future changes in the\nupwelling intensity (Fig. R5). Consistent with a recent study (Jing et al. 2023), in most\nEBUSs, the secular changes in upwelling intensity is primarily attributed to the\ngeostrophic flow changes. However, as latitudes increase, the dominant influence of\ngeostrophic flow gradually diminishes, yielding to the prominence of wind-induced\nprocesses. Specifically, in the poleward region of Canary Current system and Humboldt\nCurrent system, the wind-derived Ekman transport exceeds the geostrophic flow,\ndetermining the future changes in upwelling intensity. The intensified wind stress at\nthese high-latitude EBUSs can be attributed to the poleward migration of atmospheric\nsystems due to greenhouse warming (Rykaczewski et al., 2015). However, the future\nchanges in cross-shore geostrophic flow may be influenced by various factors,\nincluding the latitudinal variation in surface warming, asymmetrical precipitation\npatterns, and sea ice melting in polar regions (Cheng et al., 2022; Notz & Stroeve, 2016).\nConsequently, it is anticipated that distinct patterns of the cross-shore geostrophic flow\nwill manifest in individual EBUS, thereby warranting further investigation in future\nstudies.\nFig. R5. Trend of upwelling intensity calculated by UI (purple), UI (blue), and UI (red)\ng e\nderived from CESM-H. The error bar represents the standard error of linear regression.\nWe have the above analysis and discussion in the revised manuscript. See lines 156-\n175 and 233-243.\nCheng, L., von Schuckmann, K., Abraham, J. P., Trenberth, K. E., Mann, M. E., Zanna,\nL., et al. (2022). Past and future ocean warming. Nature Reviews Earth &\nEnvironment, 3(11), 776–794. https://doi.org/10.1038/s43017-022-00345-1\nJing, Z., Wang, S., Wu, L., Wang, H., Zhou, S., Sun, B., et al. (2023). Geostrophic flows\ncontrol future changes of oceanic eastern boundary upwelling. Nature Climate\nChange, 1–7. https://doi.org/10.1038/s41558-022-01588-y\nNotz, D., & Stroeve, J. (2016). Observed Arctic sea-ice loss directly follows\nanthropogenic CO2 emission. Science, 354(6313), 747–750.\nhttps://doi.org/10.1126/science.aag2345\nRykaczewski, R. R., Dunne, J. P., Sydeman, W. J., García-Reyes, M., Black, B. A., &\nBograd, S. J. (2015). Poleward displacement of coastal upwelling-favorable winds\nin the ocean’s eastern boundary currents through the 21st century. Geophysical\nResearch Letters, 42(15), 6424–6431. https://doi.org/10.1002/2015GL064694\nSong, F., Leung, L. R., Lu, J., & Dong, L. (2018). Seasonally dependent responses of\nsubtropical highs and tropical rainfall to anthropogenic warming. Nature Climate\nChange, 8(9), 787–792. https://doi.org/10.1038/s41558-018-0244-4\n4) Refinement of methodology concerning shifting definition in monthly data. A critical\naspect requiring attention in the methodology is the method of shifting definition\nemployed in analyzing upwelling duration based on monthly data, utilizing the epoch\ndifference technique (Brady et al., 2017). While this method proves effective for\nstudying the emergence of anthropogenic impacts on upwelling intensity, its\napplicability to upwelling duration analysis is questionable. As shown in Figure 2a-c,\nthe onset and termination time suffer strong interannual/decadal variabilities. From\n1920 to 2010, onset time varies from 01-15 to 03-01 in northern CalCS and Iberian\nsystem, and from 07-15 to 08-15 in Southern HCS system. Such variability suggests\nthat changes in onset time using a 30-year epoch difference can be influenced by\ninternal variability rather than anthropogenic forcing, especially when dealing with\nmonthly data. The same problem arises from the changes in termination time. This\nimplies that the comparison between CMIP6 and CESM-H, based on this methodology,\nmay not be a fair one. To address this issue and ensure a more robust comparison, it\nis advisable to focus on confirming the role of surface wind in influencing changes in\nupwelling duration. Subsequently, employing daily wind data to estimate duration in\nCMIP6 models, as demonstrated by Wang et al. (2015), can provide a more accurate\nrepresentation of the temporal dynamics of upwelling.\nThanks for your suggestion, we opt for utilizing monthly modeled vertical velocities\ninstead of the daily wind-derived upwelling index due to the following reasons. First, a\nkey motivation of this study lies in the assessment of shifting rates of upwelling onset\nand termination under greenhouse warming. However, significant disparities exist in\nestimating the onset and terminational time of upwelling season between the modeled\nvertical velocity and wind-derived upwelling index (Fig R6), which may introduce\nsubstantial uncertainties when projecting its future changes. Second, the utilization of\nthe wind-derived upwelling index for projecting future changes in upwelling relies on\nthe underlying assumption that variations in coastal upwelling are predominantly driven\nby wind stress. This assumption could hold when the time scale is shorter, as proved in\nobservations (Chereskin, 1995). However, over longer time scales, oceanic and\natmospheric processes beyond wind stress may exert a significant influence on the\nupwelling variation (Ding et al., 2023; Jing et al., 2023; Rykaczewski et al., 2015). In\nparticular, recent studies have demonstrated that oceanic processes may surpass the\ninfluence of wind stress, playing a pivotal role in modulating long-term changes in\nupwelling under greenhouse warming (Jing et al., 2023). Indeed, we have found that\nthe geostrophic transport serves as the primary driver for the seasonal advancement and\nprolonged duration of upwelling in the Atlantic EBUSs (Fig R3). Above all, it may not\nbe appropriate to solely consider wind changes when projecting the future changes in\nthe timing of the upwelling season.\nFig. R6 Climatology UI derived from vertical velocity (black) and wind (red).\nHowever, in accordance with your concerns about the influence from\ninterannual/decadal variabilities, we extend the 30-year period to a 50-year period and\nredraw Fig. 2e-h (Fig. R7). The results continue to align with the that derived from the\ndaily output in CESM-H, indicating that the observed changes could primarily represent\nthe impact of global warming.\nFig. R7 Same as Fig. 2e-h, but for 30-year (a-d) 50-year (e-h) epoch differences.\nChereskin, T. K. (1995). Direct evidence for an Ekman balance in the California Current.\nJournal of Geophysical Research: Oceans, 100(C9), 18261–18269.\nhttps://doi.org/10.1029/95JC02182\nDing, H., Alexander, M. A., & Ting, M. (2023). Revisiting the Relationship between\nthe North Pacific High and Upwelling Winds along the West Coast of North\nAmerica in the Present and Future Climate. Journal of Climate, 36(23), 8211–\n8224. https://doi.org/10.1175/JCLI-D-23-0238.1\nJing, Z., Wang, S., Wu, L., Wang, H., Zhou, S., Sun, B., et al. (2023). Geostrophic flows\ncontrol future changes of oceanic eastern boundary upwelling. Nature Climate\nChange, 1–7. https://doi.org/10.1038/s41558-022-01588-y\nRykaczewski, R. R., Dunne, J. P., Sydeman, W. J., García-Reyes, M., Black, B. A., &\nBograd, S. J. (2015). Poleward displacement of coastal upwelling-favorable winds\nin the ocean’s eastern boundary currents through the 21st century. Geophysical\nResearch Letters, 42(15), 6424–6431. https://doi.org/10.1002/2015GL064694\nMinor concerns:\n1) Descriptions of upwelling regions are incomplete. Firstly, the authors did not\nexplicitly define Northern CalCS, Iberian, and Southern HCS regions, potentially\nconfusing readers. Secondly, the authors should describe methods to calculate the trend\nerrors of onset and termination time that are shown in Figure 2a-c. Thirdly, the term\n\"equatorward regions\" is used without a precise definition.\nThanks for you comments, the upwelling regions’ definition has been added in line 433-\n443.\nThe methods we used to calculate the linear trend and their standard errors is the\nCochrane-Orcutt method (Cochrane & Orcutt, 1949), detailed description has been\nadded in line 488-496\n“The linear trends and their standard errors were obtained by regressing the time series\nof each upwelling metric against time using the Cochrane-Orcutt (C-O) method\n(Cochrane & Orcutt, 1949). This method is employed to address the autocorrelation\nwithin the time series. The model for time series of a given quantity 𝜃(𝑡) is expressed\nas 𝜃(𝑡) = 𝛽 +𝛽 𝑡+𝜖(𝑡), where the intercept (𝛽 ) and linear trend (𝛽 ) can be easily\n0 1 0 1\ncomputed based on the ordinary least-squares method. Here, 𝜖(𝑡) represents a\nstationary stochastic process with a zero mean. The standard error of the linear trend\nwas computed using the C-O method, assuming that the autocorrelation in 𝜖(𝑡) can be\napproximated as a first-order autoregressive process.”\nWe have reselected the subregions to poleward, central, and equatorward parts, the\ndefinition of the subdivision can be seen in line 422-431.\nCochrane, D., & Orcutt, G. H. (1949). Application of least squares regression to\nrelationships containing auto-correlated error terms. Journal of the American\nStatistical Association, 44(245), 32–61.\nhttps://doi.org/10.1080/01621459.1949.10483290\n2) What are the x-axis and numbers -3, -2, ….2, 3 in Figure 2d-f? Authors should\ndescribe it in the caption.\nThanks for your suggestion, the figure capture has been revised more detailed.\n“The positive values (+1, +2, and +3; warm color) indicate the months in which the\nonset and termination of the upwelling season have advanced, whereas the negative\nvalues (-1, -2, and -3; cool color) represent the months in which the onset and\ntermination of the upwelling season have delayed. The +0 and -0 (hatched) also indicate\nthe advancement and delay of the upwelling onset and termination but with the shifting\nrate less than one month. As for the duration, the positive and negative values\ncorrespond to the months that the upwelling season has prolonged and shortened,\nrespectively.”\n3) The definition of upwelling intensity in the study is based on the time and spatial\nintegral of vertical velocity at 50 m, with onset and termination times determined using\nmonthly data. However, an important aspect requiring clarification is the attribution of\nchanges in upwelling intensity. The study does not explicitly specify whether alterations\nin upwelling intensity are primarily driven by variations in upwelling duration, vertical\nvelocity at 50 m, or a combination of both factors.\nThanks for your suggestion. Following this comment, we have divided the upwelling\nintensity into two parts, the velocity-induced (the magnitude of vertical velocity at 50\nm) changes and the duration-induced changes. To ascertain the impact of velocity\nchanges on variations in upwelling intensity, we recalculated the upwelling intensity by\nkeeping the onset date and termination date fixed at their climatological values\n(represented as the velocity-induced upwelling intensity). Then, the differences\nbetween the upwelling intensity and the velocity-induced intensity indicates the\ncontribution of duration changes (duration-induced upwelling intensity).\nFor the annual upwelling regions, where the duration remains constant at 365 days,\nfuture changes in upwelling intensity are solely attributed to the changes in the velocity-\ninduced upwelling intensity. In subregions where changes in duration occur, the\nupwelling intensity is also primarily attributed to the velocity-induced upwelling\nintensity changes across all EBUSs (Fig. R10). Thus, alterations in upwelling intensity\nare primarily driven by the changes in the magnitude of the vertical velocity.\nA brief discussion has been included in the manuscript, please refer to line 232-233.\nFig. R10. Trend of changing percentage in upwelling intensity, UI-induced intensity\nand duration-induced intensity in CalCS-P, CanCS-P, CanCS-E, and HCS-P regions\nduring 1920-2100. The error bar denotes the standard error of linear regression.\n4) The authors consistently express the agreement percentages between CMIP6\nsimulations and CESM-H, but it is crucial to note that not all CMIP6 simulations were\nutilized in the paper. To enhance precision and clarity in reporting, it is recommended\nto present the results in a format that specifies the number of simulations in agreement,\nsuch as \"xx out of xx simulations.\" This format provides a more accurate representation\nof the agreement statistics, ensuring that readers have a clear understanding of the scope\nand basis for the reported percentages.\nThanks for your helpful suggestions, this has been revised in the manuscript. Specific\nnumber of simulations were added in the brackets after the use of percentage in the\nmanuscript.\n5) CESM-H models used in this study do not include biogeochemistry components,\nwhich means that nitrate in CESM-H is a passive tracer. Given this, it is not surprising\nthat changes in vertical velocity dominate changes in vertical nitrate fluxes, as no\nadditional sources and sinks of nitrate are considered in the future. To deepen the\nunderstanding of the biogeochemical dynamics, it is recommended to explore whether,\nin CMIP6 models with biogeochemistry components, changes in vertical nitrate flux\nare similarly dominated by changes in vertical velocity. This comparison will shed light\non the role of biogeochemistry in influencing vertical nitrate flux patterns. In addition,\ntotal vertical nitrate flux calculated using monthly data does not include the component\ninduced by submonthly variabdilities (equivalent to eddy component). The eddy flux is\nresolved by CESM-H but parameterized in 1° or coarser CMIP6 models. Therefore, it\nis advisable to quantify the eddy fluxes, leveraging the advantages of CESM-H, to\nprovide a more comprehensive understanding of the nuances in vertical nitrate fluxes\nand their implications.\nSorry about the misleading presentation. The vertical nutrient transport (VNT) and the\nfollowing NPP work are all evaluated using the CMIP6 ESMs, due to the lack of\nbiogeochemical components in CESM-H. We’ve added a sentence in line 249-251 for\nbetter expression.\n“To test this hypothesis, we quantify the long-term changes in the upward nutrient\ntransport and NPP based on the Earth System Models (ESMs) in CMIP6.”\nWe fully agree that leveraging CESM-H for quantifying eddy fluxes is indeed\nadvantageous. However, integrating such a high-resolution model with biogeochemical\ncomponents remains challenging. So, diagnosing the roles played by eddies that\nmodulate VNT is still an ongoing work. We added some discussion about the\nshortcomings of the CMIP6 models to caution the readers, mainly about its spatial and\ntemporal resolution, please see in line 303-306 and 311-314\nSpecific comments:\nL34: EBUSs generally stands for eastern boundary upwelling systems, not coastal\nupwelling systems.\nRevised.\nL41-42: This statement is not accurate. The reason for using wind to estimate upwelling\nis because upwelling is mainly wind-driven, not because sparsity of oceanic\nobservations.\nThanks for your comment. We have revised this sentence as\n“Conventionally, wind-induced Ekman transport has been regarded as the primary\nmechanism driving upwelling in the EBUSs (Bakun, 1973).”\nL47-48: Previous studies showed that projected changes in wind and upwelling both\nshow inconsistency across different models. But in one certain model, if it is not wind\ndetermining future changes in upwelling intensity, what are the other factors? Please\nclarify it here.\nFollowing your comment, the other factors has been introduced more specifically.\nPlease see in line 53-55.\n“In particular, future changes in stratification (Oyarzún & Brierley, 2019; Sousa et al.,\n2020; Vázquez et al., 2023) and geostrophic transport (Ding et al., 2021; Jing et al.,\n2023) could surpass the impact of winds, …”\nFigure 1a: vertical velocity at which depth?\nThanks for your suggestion, the vertical velocity at 50 m, this has been added in the\nfigure caption.\nL76: ‘Consistent with the existing theoretical arguments and simulations’, please add\nreferences here.\nThanks for your suggestion, the references were added.\nL83: The maximum between 14°S-17°S in HCS and 27°S-28°S in BCS occurs in the\nearly austral spring. Please precisely define poleward and equatorward regions for HCS\nand BCS.\nThanks for your comments, this has been revised. Please see in line 422-443 for the\nregion definition, and the revised seasonal description can be seen in line 87-97\nL86-87: Because duration of upwelling is the major focus of this paper, it will be more\ncomprehensive to provide a figure including duration as a function of latitude in each\nupwelling system. It is difficult to follow the authors’ argument without any figures like\nthat.\nThanks for your suggestion. We have displayed a figure (Figure R11) imitating the\nFigure 1 in Wang et al. (2015) to illustrating the latitudinal distribution of duration. This\nhas been added to the Supplementary Fig. 1.\nFig. R11. Latitudinal distribution of upwelling duration during 1920-1949. The shading\ndemonstrating the standard error across the 30 years coverage.\nL104: ‘all EBUSs’ means Northern CalCS, Iberian, and Southern HCS regions or\nregions defined by boxes in Figure 1a?\nSorry for the misleading description, this has been avoided in the revised manuscript.\nL129: Authors divided EBUSs into poleward and equatorward regions, which show\ndifferent changes in duration (Figure 2 and Supplementary Figure 2). ‘The consistent\nchanges across all EBUSs …’ is confusing.\nSorry for the misleading description, we have revised the word “consistent” and “all”\nin the manuscript.\nL123-124: 75% and 73% are equivalent to how much out of how much?\nRevised.\nL133: ‘North Atlantic’?\nSorry about the mistake, this has been revised.\nL134-135: Canary and Iberian regions are not defined.\nRevised, please see in line 432-443\nL138: north and south portions of CalCS are not defined.\nRevised, please see in line 432-443\nL140-141: As shown in Figure 1, the Benguela system extends from 33°S to 15°S.\nAuthors cannot simply group it into equatorward systems.\nThanks for your comments, we have reevaluated the subdivision of the BCS, and the\ndescription of BCS upwelling changes has been revised accordingly.\nL157, 159, 161: Please add actual numbers (xx out of xx) after the percentage because\nyou are not using all models that participated in CMIP6.\nAdded.\nL336: Should be ‘Brady et al.’\nSorry about the misspelling, this has been revised.\nReferences:\nWang, D., Gouhier, T. C., Menge, B. A., & Ganguly, A. R. (2015). Intensification and\nspatial homogenization of coastal upwelling under climate change. Nature, 518(7539),\n390-394.\nBrady, R. X., Alexander, M. A., Lovenduski, N. S., & Rykaczewski, R. R. (2017).\nEmergent anthropogenic trends in California Current upwelling. Geophysical Research\nLetters, 44(10), 5044-5052.\nRykaczewski, R. R., & Dunne, J. P. (2010). Enhanced nutrient supply to the California\nCurrent Ecosystem with global warming and increased stratification in an earth system\nmodel. Geophysical Research Letters, 37(21).\nREVIEWER COMMENTS\nReviewer #1 (Remarks to the Author):\nThe authors have done an excellent job addressing each of the questions raised in the review and\nimplementing the various suggestions made by the reviewers. After reading the responses, I don't\nhave any further doubts to resolve. Therefore, I propose that the article be accepted for publication\nin the journal.\nReviewer #3 (Remarks to the Author):\nThis is my second review of this paper. I appreciate the authors’ efforts to enhance the introduction,\nrefine the definition of upwelling regions, and integrate the discussions on the contributions of\nwinds and geostrophic flow to upwelling changes. However, I cannot recommend the publication of\nthis paper in Nature Communications with the current version. The reasons for this decision are\nstated as follows.\n1. Some results are contradictory to previous studies. Authors show a decreasing trend in upwelling\nintensity in the central BCS. However, previous studies (Wang et al., 2015; Chang et al., 2023)\ndemonstrated an increasing trend.\n2. Some evidence in the paper does not support the corresponding arguments.\n(1) Authors argue that results within 50 km from the coast are similar to these within 200 km. This is\nnot entirely true. Firstly, the detected termination and onset dates are totally different, with the time\ndiscrepancy of about one month between these two cases. This discrepancy is important for\npolicy-making in the upwelling regions. In addition, the trend of upwelling intensity in BCS-C is -\n1.6±0.7% per century using 50 km, but -7.1±0.6% per century using 200 km. The former one is hard\nto identify as a significant trend compared to the latter one.\n(2) I am still not convinced that the CMIP6 results can support CESM-H results in terms of the\nchanges upwelling durations (Fig. 2 Monthly data granularity cannot adequately capture variations\nof less than one month. However, using monthly data, authors identified that the changes in\nupwelling termination and onset dates are less than one month in some CMIP6 models. These\nresults are questionable. Additionally, the changes in upwelling duration in CESM-H are ~10.8, 21.2,\n38.3, 19.3 days per century in CalCS-P, CanCS-P, CanCS-E, and HCS-P, respectively (Fig. 2). These\nchanges are mostly less than one month, which is impossible to be supported by monthly data.\n3. The story is incomplete. The upwelling biases in CEMS-H are discussed (Fig. 1) in the beginning of\nthe paper, lending confidence to the future projections in the upwelling duration and intensity (Fig. 2\nand Fig. 3). To be a complete story, the discussion of impacts of upwelling in productivity should\nalso use CESM-H, instead of only CMIP6 models because upwelling biases in CMIP6 are not\ndiscussed in this paper. If there is no biogeochemistry component in CEMS-H, I do not think it is\nappropriate to include the productivity results as a primary section in this work because readers\nwill question the reliability of these findings.\n4. The definitions of land-sea temperature (pressure) differences and atmospheric pressure in the\ninterior ocean are not defined in Supplementary Figure 6. Additionally, x-axis in Figure 2e-h is still\nnot introduced.\nReply to the third reviewer\nWe are very grateful to you for your time in carefully reading our revised\nmanuscript and providing helpful comments that make our manuscript better. We have\ncarefully considered each of your comments (in blue) and revised the manuscript\naccordingly. Please find our response (in black) to your comments below.\nThis is my second review of this paper. I appreciate the authors’ efforts to enhance\nthe introduction, refine the definition of upwelling regions, and integrate the discussions\non the contributions of winds and geostrophic flow to upwelling changes. However, I\ncannot recommend the publication of this paper in Nature Communications with the\ncurrent version. The reasons for this decision are stated as follows.\n1. Some results are contradictory to previous studies. Authors show a decreasing trend\nin upwelling intensity in the central BCS. However, previous studies (Wang et al., 2015;\nChang et al., 2023) demonstrated an increasing trend.\nThanks for your comments. As you pointed out, the upwelling intensity in the central\nBCS exhibits a decreasing trend, which seems to contradict previous findings by Wang\net al. (2015) and Chang et al. (2023) (hereafter referred to as W15 and C23,\nrespectively). Upon thorough investigation, we have identified the reasons for these\ndifferences in upwelling intensity trends.\nThe differences between our study and W15 arise primarily from variations in the data\nutilized to calculate upwelling intensity. W15 determined upwelling intensity using the\nwind stress, whereas our study used model output vertical velocity. Consequently, W15\nattributed secular variations in upwelling intensity solely to Ekman transport, whereas\nour study considers changes from both the Ekman transport and geostrophic transport.\nIn order to directly compare our findings with W15, we further decompose the changes\nin upwelling intensity into components associated with variations in wind and\ngeostrophic flow. In particular, the wind-induced upwelling intensity (UI , red bar) in\ne\nBCS-C in our study shows an increasing trend, consistent with the findings in W15.\nHowever, the increasing trend is counteracted by a decrease in geostrophic-induced\nupwelling intensity (UI , blue bar), resulting in an overall negative upwelling intensity\ng\ntrend in BCS-C. Therefore, our findings are not contradictory to those of W15. Indeed,\nour study suggests that the changes in geostrophic flow may exert a significant\ninfluence on future upwelling intensity. We have added above discussion and highlight\nthe import role of geostrophic transport on future upwelling intensity changes in the\nrevised manuscript. See lines 241-244.\nFig. R1 Upwelling intensity percentage trends calculated by vertical velocity (purple),\nwind-induced upwelling index (UI , red) and geostrophic-induced upwelling index (UI ,\ne g\nblue).\nIn contrast, the disparity in future upwelling intensity changes between our study and\nC23 primarily stem from the variation in offshore distance chosen for averaging. In\nparticular, C23 focused on future changes in coastal upwelling within a 50-km-wide\ncoastal zone, while we employed a 200-km average that includes both coastal upwelling\nand a large fraction of offshore wind stress curl-driven upwelling. Following C23,\nwe\nhave recomputed the upwelling intensity for the approximately 50-km-wide coastal\nzone in the BCS regions. The secular changes in both poleward and equatorward\nregions exhibit a resemblance to those calculated at a distance of 200 km from the\ncoastline (Fig. R2a and c). However, notable differences are evident in the central\nregion. Specifically, we observed a negative trend in upwelling within the 200-km-wide\nzone across nearly all latitudes, but a distinct positive trend at specific latitudes within\nthe 50-km-wide coastal zone, similar to the findings in Figure 3j of C23. This indicates\nthat the secular changes in upwelling intensity exhibit a divergent pattern at varying\ndistances from the shore in the BCS-C region.\nThe contour plot (Fig. R2b) provides a clearer depiction of the offshore distribution\nwithin BCS-C, showing the most significant decreases in upwelling intensity lays\nbetween 50-100 km offshore, with weaker positive trends between 20-50 km offshore.\nTherefore, the varying cross-shore distribution in this region is likely to be the primary\nfactor contributing to the contrasting results between our study and C23. However, it is\nworth noting that the secular changes in upwelling intensity within the 50-km-wide\ncoastal zone in most EBUSs exhibit similarities to those calculated at a distance of 200\nkm from the coastline. Nevertheless, we also caution readers about considering the\ncross-shore variations in trends in upwelling intensity in some certain EBUSs,\nparticularly in BCS-C region. See lines 199-211.\nFig. R2 (a) Latitudinal distribution of upwelling intensity trend averaged within 200 km\noffshore. (b) Distribution of upwelling intensity varies with latitude and offshore\ndistance. (c) Same as (a) but of 50-km average.\n2. Some evidence in the paper does not support the corresponding arguments.\n(1) Authors argue that results within 50 km from the coast are similar to these within\n200 km. This is not entirely true. Firstly, the detected termination and onset dates are\ntotally different, with the time discrepancy of about one month between these two cases.\nThis discrepancy is important for policy-making in the upwelling regions. In addition,\nthe trend of upwelling intensity in BCS-C is -1.6±0.7% per century using 50 km, but -\n7.1±0.6% per century using 200 km. The former one is hard to identify as a significant\ntrend compared to the latter one.\nThanks for your comments. We fully agree with you that there is a discrepancy in the\ndetected onset and termination dates between regions located within a 50-km-wide band\nfrom the coastline and those situated within a 200-km-wide band. Indeed, the\nclimatological onset and termination date of the upwelling season changes with distance\nfrom the coast (Figure R3). However, their secular changes under greenhouse warming\nexhibit consistent values, thereby providing robust support for the validity of our\nmanuscript's conclusions.\nRegarding the intensity of upwelling, the long-term variations in upwelling intensity\nwithin the 50-km-wide coastal zone in most EBUSs exhibit similarities to those\ncalculated at a distance of 200 km from the coastline, albeit with distinct magnitudes.\nHowever, we also acknowledge that the trends of upwelling intensity could exhibit a\ndivergent pattern across cross-shore direction, particularly in BCS-C region (see our\nreply to your comment 1).\nIn this study, we seek a first-order understanding of upwelling timing and intensity\nchange in EBUSs under global warming, which has not been previously addressed in\nrecent studies. In addition, integrating vertical velocity over a sufficiently wide region\nshould be less sensitive to model resolution and more suitable for the CMIP6 analysis.\nFor these reasons, we focus on the upwelling changes averages in a 200-km-wide band\nfrom the coast. Such integration covers coastal upwelling and a large fraction of\noffshore wind stress curl-driven upwelling, both of which are suggested to influence\nthe ecosystem. We know that the comprehensive understanding of the detailed structure\nof upwelling changes within individual EBUSs is undoubtedly important, but it exceeds\nthe scope of this study and warrants further investigation in future research.\nWe have added Figure R3 and briefly discussed the results in the revised manuscript.\nPlease see Supplementary Figure 4, lines 130-134 and 199-211.\nFig. R3 a-d, The climatological and future upwelling onset and termination date in\nCalCS-P, CanCS-P, CanCS-E and HCS-P. The x-axis denotes the average distance from\nthe shore.\n(2) I am still not convinced that the CMIP6 results can support CESM-H results in terms\nof the changes upwelling durations (Fig. 2 Monthly data granularity cannot adequately\ncapture variations of less than one month. However, using monthly data, authors\nidentified that the changes in upwelling termination and onset dates are less than one\nmonth in some CMIP6 models. These results are questionable. Additionally, the\nchanges in upwelling duration in CESM-H are ~10.8, 21.2, 38.3, 19.3 days per century\nin CalCS-P, CanCS-P, CanCS-E, and HCS-P, respectively (Fig. 2). These changes are\nmostly less than one month, which is impossible to be supported by monthly data.\nThanks for your comments. We admit that monthly data may not sufficiently capture\nprecise variations in the onset and termination of upwelling, particularly for changes\noccurring within a time frame shorter than one month. However, due to storage\nlimitations, daily output of vertical velocity is not provided in all CMIP6 climate\nmodels. In this study, our objective is to provide a precise characterization of the shift\nrates in the onset and termination of upwelling seasons using daily output from the\nCESM-H model. Additionally, analysis from CMIP6 serves as supplementary evidence\nto support the primary conclusion derived from the CESM-H model. Indeed, numerous\nstudies have utilized monthly data to investigate the seasonal timing changes in\nvariables such as SST (Allan & Allan, 2019; Dwyer et al., 2014; Liu et al., 2024),\nprecipitation (Feng et al., 2013; Song et al., 2018, 2021), and also upwelling (Brady et\nal., 2017; Oyarzún & Brierley, 2019; Rykaczewski et al., 2015).\nIn this manuscript, we assessed the changes of upwelling onset and termination by\nexamining the variations in upwelling during the transitional month, which has been\nwidely used in previous studies (Brady et al., 2017; Oyarzún & Brierley, 2019;\nRykaczewski et al., 2015). To demonstrate that monthly data can reflect changes in\nupwelling timing, we conducted further analysis using three other commonly used\nmethods:\n1. Linear interpolation: linearly interpolate the monthly data to generate daily data and\nassess the onset and termination date of upwelling season.\n2. Polynomial fitting: fit a 4th-order polynomial to the monthly data and use the\npolynomial to obtain the onset and termination date of upwelling season.\n3. First harmonics: performing Fourier analysis on the monthly data and extracting the\nfirst harmonic mode. It is worth noting that the mean value is also important for the\nupwelling timing detection, therefore, not only the first harmonic mode but also the\nmean value (harmonic mode with zero frequency) was also included in this method\nto accurately determine upwelling timing.\nAs shown in Fig. R4, despite slight deviations observed in the reconstruction of precise\nupwelling onset and termination dates by each method, all approaches qualitatively\nreflect the secular changes derived from the model's daily output, with discrepancies\nless than 3 day per century for both onset and termination in most regions. We further\napply these methods to the monthly output from the CMIP6 models (Fig. R5). Most\nCMIP6 models exhibit a seasonal advancement and extended duration of upwelling,\nand showing consistency across different reconstruction methods. This leads further\nsupport to our conjecture.\nFig. R4 a-c, Examples of the climatological results derived from direct daily output\n(black line), monthly output (blue dots), and reconstructed daily results via different\nmethods (colored lines). d-g, Projected changes in upwelling timing derived by daily\nCESM-H output (black) and monthly CESM-H outputs using various methods (colors)\nin CalCS-P, CanCS-P, CanCS-E, and HCS-P.\nFig. R5 Projected changes in upwelling onset date (a,d,g), duration (b,e,h), and\ntermination date (c,f,i) derived from monthly CMIP6 models in CalCS-P, CanCS-P,\nCanCS-E, and HCS-P. Each row represents the results from different methods.\n3. The story is incomplete. The upwelling biases in CEMS-H are discussed (Fig. 1) in\nthe beginning of the paper, lending confidence to the future projections in the upwelling\nduration and intensity (Fig. 2 and Fig. 3). To be a complete story, the discussion of\nimpacts of upwelling in productivity should also use CESM-H, instead of only CMIP6\nmodels because upwelling biases in CMIP6 are not discussed in this paper. If there is\nno biogeochemistry component in CEMS-H, I do not think it is appropriate to include\nthe productivity results as a primary section in this work because readers will question\nthe reliability of these findings.\nThanks for your suggestion. As you pointed out, the productivity in the EBUSs is\nheavily reliant on simulated upwelling, which may introduce significant uncertainties\nin low-resolution CMIP6 models due to their inherent biases in upwelling simulation.\nHowever, as highlighted by other reviewers, comprehending the response of ecosystem\nproductivity in the EBUSs to future climate changes is a pivotal inquiry. In this study,\nwe have unveiled an unexplored facet in the association between forthcoming\nalterations in upwelling and primary productivity under global warming. For this reason,\nwe intend to retain this part in the manuscript. However, following your comments, we\nhave reorganized the paper by moving the discussion on the implications of upwelling\nchanges on primary productivity from the main text to the discussion section, please\nsee in line 263-292.\nAdditionally, a brief discussion on the uncertainties associated with utilizing low-\nresolution CMIP6 data are included in line 303-311.\n4. The definitions of land-sea temperature (pressure) differences and atmospheric\npressure in the interior ocean are not defined in Supplementary Figure 6. Additionally,\nx-axis in Figure 2e-h is still not introduced.\nSorry for the negligence, we have addressed the methodological details in the figure\ncaption of Supplementary Fig. 6. Additionally, the x-axis for Figure 2e-h has been\nintroduced both in the figure and their captions.\nREVIEWERS' COMMENTS\nReviewer #3 (Remarks to the Author):\nThe authors have addressed all my concerns. I will recommend the publication of this work in\nNature Communications."}]

Decision letter