Introduction

Regions of strong coastal upwelling occur on the eastern edge of most ocean basins, where alongshore equatorward winds drive offshore ocean transport near the surface and replenishment with water from below1,2,3. The replenished water, higher in nutrients4, supports photosynthesis, driving thriving ecosystems and productive fishing industries. While Eastern Boundary Upwelling Systems (EBUSs) occupy only about 1% of the ocean’s surface, they produce around 20% of the global fish catch5 estimated to be worth about $400 billion by the Food and Agricultural Organization6. It is thus important to understand how EBUSs will respond to anthropogenic climate change.

Many studies have attempted to address this issue using both historical observations and IPCC climate model projections7,8,9,10,11,12,13,14,15,16. One well-known proposed mechanism, the so-called Bakun hypothesis17, postulates that under future warming the land masses will warm more than the oceans, while the increased land-ocean thermal contrast will produce stronger pressure gradients, alongshore winds, and thus more upwelling. A meta-analysis of published literature on observational data found that most EBUSs have experienced increasing wind trends in summer, especially at the higher latitudes, but not necessarily in the annual mean7. However, it is also worth noting that the winds used in Bakun hypothesis17 were biased due to a change in the size (height) of merchant ships in 1960s18, and a more recent study19 using corrected winds in the Peru Upwelling system found no upwelling favorable wind trend. Additionally, short record lengths and natural climate variability can further add large uncertainties to the analysis of future upwelling changes from historical data.

Climate models offer an alternative approach to study this problem. One analysis of alongshore wind-stress from Coupled Model Intercomparison Project phase 5 (CMIP5) models showed responses similar to Bakun’s hypothesis for the Peru-Chile (P-CUS), Canary (CUS), and Benguela Upwelling Systems (BUS) but less clearly for the California Current System10 (CCS). Another interpretation of the model outputs was that upwelling strengthened on the poleward edge, possibly due to poleward expansion of subtropical cells8. In contrast, a recent study11, using a large ensemble of the Community Earth System Model (CESM) to separate out forced signals from natural variability for the CCS region, found that upwelling intensified only in the Northern CCS in spring, whereas in summer it actually decreased. Furthermore, the study showed that anthropogenic climate change signals only emerged in the second half of the 21st century, a result that challenges trend analyses based on historical records.

However, a major issue with interpretating CMIP model projections of future upwelling change is that these models have consistently shown significant warm SST biases in EBUSs20,21,22,23, indicating poorly represented key upwelling processes. Among these are the observed atmospheric low-level coastal jets (LLCJs) forming at the top of the planetary boundary layer along the upwelling coasts with a core that is located within ~0.5° of the coast24,25,26,27,28. Typical CMIP models, with ~1° horizontal resolution, simulate a weak jet core several degrees away from the coast28,29,30,31,32. The failure of coarse resolution models to simulate the realistic jet core and the rapid drop-off of the winds towards the coast produces an erroneous representation of the wind stress curl, resulting in an exaggerated poleward heat transport from the tropics, and an underrepresented upwelling-driven vertical cooling, leading to the warm SST biases32,33,34. Meanwhile the fine-scale ocean fronts, eddies, filaments and near-coast upwelling35 are also not explicitly represented by CMIP models, while their parameterization shows large uncertainties36. All these issues call into question the suitability and creditability of standard CMIP models for projecting future changes of EBUSs. Although dynamic downscaling with high-resolution regional models can help alleviate some issues, many existing downscaling efforts use regional ocean-only models14,15,16,37 and cannot address the near-coast wind bias. Some studies have attempted to reduce the wind bias using statistically downscaled wind products14,15,16, but this modeling approach can also suffer from dynamical inconsistency and errors arising from using low-resolution global model simulations as boundary conditions.

Results

Sea surface temperature in EBUS

To address this and other societally relevant issues, we conducted an unprecedented ensemble of high-resolution historical and future climate simulations (hereafter HR) using CESM with a nominal horizontal resolution of 0.25° for the atmosphere and land components and 0.1° for the ocean and sea-ice components. The results are compared to a parallel set of low-resolution simulations (hereafter LR) at a nominal resolution of 1° for all model components38 (also see Methods). Improvements in HR over LR are evident. HR realistically simulates the narrow structure and intensity of coastal upwelling in all EBUSs, including CUS and CCS in the Northern Hemisphere (NH), and BUS and P-CUS in the Southern Hemisphere (SH) (Fig. 1), based on a comparison between  observed and simulated annual mean sea-surface temperature (SST) in these regions during 1991-2020. Although the narrow upwelling zone along the northern Benguela coast (north of ~25°S) in the BUS is still underrepresented by HR, it is far more realistic than in LR, especially in the central Benguela, and in the southern Benguela (south of ~30°S) HR even shows cooler-than-observed SST that suggests an overestimate of the upwelling. In all other EBUSs the HR-simulated SSTs show a remarkable agreement with the observations, while the narrow intense coastal upwelling is virtually absent in LR. These results hold not only for the CESM simulations, but also for other climate model simulations within the High-Resolution Model Intercomparison Project39 (HighResMIP) (Supplementary Fig. 1). Consistent with previous studies32,34, the vastly improved representation of the coastal upwelling in HR is closely linked to the better resolved LLCJs and associated wind stress curls along the upwelling coasts, particularly in terms of the jet core structures and their proximity to the coast (Supplementary Fig. 2). Since the atmospheric model resolution in HighResMIP ranges from 0.25° to 0.5° (Supplementary Table 1), it suggests that a 0.5° or finer atmospheric grid is required to simulate low-level coastal jets, and the associated coastal wind stress and wind stress curl32,40. In addition, ocean eddies and upwelling fronts are explicitly represented in HR rather than parameterized as in LR. In fact, a comparison between satellite observed and HR-simulated eddy kinetic energy (EKE) in EBUSs shows an overall agreement, except that in the Canary and Chilean region HR tends to underestimate EKE (Supplementary Fig. 3). Together, these results suggest that HR is more credible and represents dynamical processes much more realistically for projecting future changes along EBUSs. In the following, we focus our analyses on the projected EBUS changes in HR.

Fig. 1: Eastern boundary upwelling systems (EBUSs) and sea-surface temperature (SST).
figure 1

Annual mean SST bias (°C, model minus observation) from HR and locations of five EBUSs (black boxes) (a) and difference between LR and HR SST (°C, b). Observed SST45 (°C) in CCS (c), CUS (d), P-CUS (e, f), BUS (g). HR SST (°C) in CCS (h), CUS (i), P-CUS (j, k), BUS (l). LR SST (°C) in CCS (m), CUS (n), P-CUS (o, p), BUS (q). Observed and simulated SSTs are all averaged for the period of 1991–2020. Contours in (a) and (b) are −3, −1 (dashed) and 1, 3 °C (solid).

Full size image

Projections of winds, ocean temperature, and vertical velocity

To start, we examine the difference in alongshore wind stress strength between future and historical periods, since it lies at the heart of the Bakun hypothesis. From Fig. 2a–e, it is immediately clear that not all regions show strengthened winds. In the NH, the CCS shows unambiguous future weakening of the wind stress with a maximum change near 38°N just south of Cape Mendocino. In contrast, the change in the CUS is more complex; although there is an overall increase in the future alongshore wind stress, it weakens marginally near the coast. Such a change in wind stress pattern will lead to a change in near-coast wind stress curl (Fig. 2g). Consistent with previous studies8, the wind stress change in the CUS shows a poleward shift. In the SH, however, the winds off southern Chile show a substantial future strengthening, while there is also a general strengthening off the Benguela coast. The Peruvian system, however, does not exhibit significant wind stress changes (Fig. 2c). Therefore, there is clearly a regional dependence as to whether the alongshore winds follow the Bakun hypothesis. In the NH EBUSs, the coastal winds show a decrease in strength, which is at odds with the Bakun hypothesis, whereas in the SH EBUSs the winds are strengthening as predicted, except off the Peruvian coast where wind changes are mostly insignificant. Consistent with the wind stress changes, the project wind stress curl changes show a substantial future decrease in the CCS and increase off Chilean and Benguela coast, but insignificant change off Peruvian coast (Fig. 2f–j).

Fig. 2: Projected future changes of alongshore wind stress and ocean temperature in each EBUS from HR.
figure 2

Projected change (defined as the mean over the future period of 2071–2100 minus the mean over the historical period of 1991–2020) of alongshore wind stress (×10−2 N m−2) shown in longitude-latitude plane (upper) for CCS (a), CUS (b), P-CUS (c, d), BUS (e), and projected change of wind stress curl (×10−7 N m−3) (middle) for CCS (f), CUS (g), P-CUS (h, i), BUS(j). The upwelling favorable wind stress is southward (negative) in the NH and northward (positive) in the SH. Upwelling favorable wind stress curl is positive in the NH and negative in the SH. Areas without green hatching represent the regions where the mean alongshore wind stress over the historical period of 1991–2020 is upwelling favorable (see Supplementary Fig. 2). In order to avoid confusion in the sign of future change, future strengthening (weakening) in the upwelling favorable alongshore wind stress (a–e) and curl (f–j) is shown as positive (negative) values here, regardless of the EBUS hemisphere. Projected change of ocean temperature (°C in color) and historical annual mean temperature (°C in contours) shown in cross-shore section as a function of distance (°) from the coast and depth (m) (bottom) for CCS (k), CUS (l), P-CUS (m, n), BUS (o). The cross-shore sections of temperature in (k)–(o) are derived by averaging the temperature along each upwelling coast within the area indicated in green in (a)–(e) and (f)–(j). Areas without gray hatching represent the regions where the changes are significant at 95% level based on a Student’s t test.

Full size image

Next, we examine alongshore variations of future coastal upwelling changes in EBUSs. Historically, coastal upwelling has been estimated using the so-called Bakun index2. A revision to the Bakun index has been recently proposed that includes additional upwelling dynamics35. Here, we include all the dynamic processes affecting coastal upwelling in the model by computing an upwelling index using the simulated vertical velocity at the Ekman depth (see Methods) within a narrow coastal zone (Fig. 3). Comparing HR and LR shows that the HR upwelling index for the historical period is a factor of 3-5 times stronger than that in LR and shows a much richer alongshore variability (Fig. 3 left). The CUS displays the strongest mean upwelling (2.4–7.9 × 10−5 m s−1) and the strongest alongshore variability off the coast of Morocco, which appears to be linked to the small-scale topography (Supplementary Fig. 4). However, such features are completely missing in LR that shows only a gradual northward decrease in upwelling strength from 0.6 × 10−5 m s−1 in the southern CUS to 0.4 × 10−5 m s−1 in the northern CUS. Similar differences between HR and LR are also seen in other EBUS. The projected future upwelling changes show equally striking differences between HR and LR (Fig. 3 right). While LR projects generally weak and nearly uniform future upwelling changes along the coast, which agrees with previous studies10, HR projections are much more dynamic and spatially variable. For example, while LR projects a decrease in upwelling of 10−6 m s−1 along the Peruvian coast, HR projects a local maximum upwelling decrease of one order magnitude larger near 15°−16°S in the southern Peruvian coast. Along the Chilean coast, LR again projects a nearly uniform future increase of coastal upwelling, but HR shows that the upwelling increase near 37°S is 2–3 times larger in other regions. In the northern BUS between the North and Central Namibian Upwelling Cells, LR projects a future upwelling decrease, while HR projects a substantial increase of 1–2 × 10−6 m s−1. The upwelling system showing the most consistency between HR and LR is the CCS, where HR and LR both project a nearly uniform future upwelling decrease south of 40°N, except that the magnitude of the projected decrease in HR is 2–3 times of that in LR. As in the mean upwelling, the CUS stands out as showing the strongest HR-projected future upwelling changes, as well as the strongest alongshore variations. The projected alongshore variation is so strong that the sign of the future upwelling changes can alter within a short distance. This suggests that CUS is a very dynamic system, with future changes likely involving complex multiscale interactions between coastal circulations, ocean eddies, and local fine-scale topography.

Fig. 3: Upwelling index and the associated futures change off the coast of each EBUS.
figure 3

Mean upwelling index (×10−5 m/s) averaged over the historical period of 1991-2020 (left) in HR (gray shade and purple with squares) and LR (green with triangles) and the corresponding future change (×10−5 m/s) computed as the future mean (over 2071–2100) minus the historical mean (over 1991–2020) (right) in HR (dots and purple line) and LR (green with triangles), respectively, as a function of latitude for CCS (a, b), CUS(c&d), P-CUS (e, f, g, h), and BUS (I, j). The upwelling index is computed based on simulated vertical velocity at the Ekman depth averaged in a coastal bin of 1° × 1° for LR (green with triangles) and ~0.5° (cross-shore) × 0.1° (alongshore) for HR (gray shade in left panels and color dots in right panels) to reflect model resolution difference. Red and blue dots indicate future increase and decrease in the HR upwelling index, respectively, with dark (light) colors indicating passing (failing) Student’s t test with 95% significance level. A smoothed upwelling index and future change in HR (purple) are also shown by averaging over 1° bin alongshore to compare directly with LR.

Full size image

While there is a dynamical consistency between projected changes in the upwelling index and alongshore wind stress and wind stress curl in most EBUSs except off the Peruvian coast where wind changes are largely insignificant (Fig. 2c, h), future ocean temperature changes show a more perplexing picture. Here, we focus on the annual-mean vertical temperature structure averaged along each upwelling coast, as the annual-mean changes are similar to upwelling-season mean changes (Supplementary Fig. 5). In all EBUSs, the simulated vertical temperature structures during the historical period show upper ocean isotherms tilted upwards towards the coast (solid contours in Fig. 2k–o), indicative of a well simulated upwelling in these regions consistent with the simulated SST (Fig. 1). Like the winds, the projected temperature changes again show differences between the NH and SH (color in Fig. 2k–o). In the NH, the CCS shows a strong and broad upper-ocean warming, but within 100 km (~1°) of the coastal zone where the upwelling is the strongest, future warming is reduced in the upper 50 m (Fig. 2k). The CUS is also projected to have less future warming near the coast compared to the open ocean, with the projected cross-shore gradient change being more prominent than in the CCS (Fig. 2l). These projected temperature changes seemingly imply strengthened Ekman upwelling, consistent with the Bakun hypothesis, although the coastal winds in these regions are projected to weaken (Fig. 2a, b), inconsistent with strengthened Ekman upwelling. In the SH, all three EBUSs show enhanced warming within the coastal zone in response to anthropogenic forcing, which clearly contradicts the Bakun hypothesis, even though the coastal wind changes (Fig. 2d, e) and the upwelling index changes (Fig. 3h–j), especially for Chile and the Benguela, are in line with Bakun’s ideas. These results hold not only for the CESM, but also for the multi-model ensemble mean of HighResMIP models (Supplementary Fig. 6).

Heat budget for the ocean temperature increase

Clearly, there is a discrepancy between temperature in the upwelling region and coastal wind responses to future warming as portrayed by the Bakun hypothesis that is entirely based on vertical Ekman upwelling changes. This discrepancy is, to a significant extent, caused by neglecting the contribution from horizontal transport of heat from the tropics to EBUS in Bakun’s mechanism. Figure 4 shows the 2000-2100 warming trend in the upper 50 m temperature, area-averaged within the ~1° coastal zone of the EBUSs, along with contributions from each term in the heat budget (see Methods). The warming trend varies from 0.16 Wm−2 in the CUS to 0.27 Wm−2 in Peru and CCS, and the cause of the warming differs drastically between NH and SH EBUSs. For all SH EBUSs, ocean heat advection is the dominant contributor to the warming trend, while surface fluxes are the dominant terms counteracting the warming. Decomposing the heat advection into its components indicates that it is the mean horizontal advection that makes the primary contribution to warming in the SH EBUSs (Fig. 4b). The warming trend from the mean horizontal advection can arise from 1) equatorial oceans, particularly the equatorial Pacific, warming faster, and 2) stronger LLCJs leading to an increase in the upwelling-favorable wind stress curl near the coast, which in turn drives stronger poleward coastal currents through Sverdrup balance. Both processes enhance warm advection from the tropics to the upwelling regions, producing faster warming within the coastal zone than in the open ocean. Because the Peruvian system is closest to the equator and shows a weak wind stress curl change (Fig. 2h), the warm advection in this region likely comes primarily from the faster warming equatorial Pacific. Indeed, HR shows a prominent future warming of SST in the equatorial Pacific and to a lesser extent warming in the equatorial Atlantic (Supplementary Fig. 7). On the other hand, the Chilean and Benguela systems see significant strengthening in the local wind stress curl under future warming (Fig. 2i, j), suggesting that the intensified alongshore poleward flows contribute more to the warm advection. Neither of these warming mechanisms are included in the Bakun hypothesis and can only be fully captured by high-resolution climate models, because LLCJs and the associated coastal wind stress curl, both poorly represented in low-resolution models, are vital. We emphasize that the key difference between HR and LR occurs near the coast. In LR, the warm advection via Sverdrup transport is too broad and very weak near the coast. Therefore, it has little impact on the coastal warming, which is sharp contrast to the strong and narrow coastal warm advection in HR.

Fig. 4: Heat budget for the warming trend in each EBUS.
figure 4

a The warming trend (red), and contributions from surface heat flux (yellow), total ocean heat advection (cyan), and turbulent mixing (blue). b Contributions from eddy-induced advection (green) and mean-current-induced advection (olive green) to total ocean heat advection (cyan in (a)), and contributions from horizontal mean advection (seagreen) and vertical mean advection (gray) to mean-current-induced advection (olive green). c Contributions from shortwave (green), net longwave (olive green), latent (seagreen), and sensible (gray) heat flux to the total heat flux (yellow in (a)). The trend in each budget term (in W m−2) is computed based on averages within the top 50 m and over a 100 km wide strip next to the coast within each green area indicated in Fig. 2a–e, using a 3-member ensemble mean of monthly mean output from HR for the period of 2000–2100 after subtracting the preindustrial control simulation. Note that there is a large amount of cancellation between the advection and surface heating terms, such that the net trend is relatively small compared to these individual terms, but remains positive for all EBUSs.

Full size image

In NH the largest contributor to the warming trend is the net surface heat flux, indicating that the coastal warming trend is driven by atmospheric heating. Ocean heat advection is primarily responsible for driving a cooling trend, in sharp contrast to the SH. The cold advection is ~3 times stronger in the CUS than in the CCS (Fig. 4). Further decomposition of the advection shows the mean vertical heat advection, rather than the horizontal advection as in the SH, is the dominant contributor to the cooling trend, and its strength in the CUS is more than twice of that in the CCS. This stronger vertical heat advection cannot be explained by increased vertical Ekman upwelling because the coastal winds are weakened under future warming (Fig. 2). As the vertical velocity in the CCS decreases (Fig. 3a, b) in the future, the enhanced cooling due to the vertical heat advection here is most likely caused by the increased vertical temperature gradient within the upwelling zone as the increased surface heat flux increases upper ocean stratification. This differs from the CUS, where despite a decrease in the alongshore wind stress, the vertical velocity shows a future increase (Fig. 3c, d). This increase can be attributed to the intensification of near-coast wind stress curl as the wind pattern changes off the CUS under future warming (Fig. 2). Indeed, there is a narrow stretch of strengthened wind stress curl along the Canary coast from 25°N to 28°N (Fig. 2g), which coincides with the large increase of the near-coast vertical velocity (Fig. 3). These results suggest that the strong increase in the vertical cold advection in the CUS depends on increases in both ocean stratification and vertical velocity driven by intensified near-coast wind stress curl through Ekman suction. Both these mechanisms are more complex than the classical Ekman upwelling mechanism invoked by Bakun.

Regarding the contrasting role of surface heat fluxes in the NH and SH, a decomposition of surface heat fluxes shows that the shortwave is always decreasing due to the aerosol scattering and absorption41 (Fig. 4c). In contrast, the net longwave radiation is seen to help warming in all EBUSs. Interestingly, turbulent heat fluxes show a different role between NH and SH, with less heat release from the ocean in NH but more heat release in SH. The changes in the turbulent heat fluxes are consistent with the wind changes, i.e., less turbulent heat loss in NH where the winds get weaker and more in SH, in agreement with the increasing winds. Thus, turbulent heat fluxes work in concert with longwave radiation to drive warming in the NH EBUSs, while in the SH EBUS they act with shortwave radiation to drive cooling .

Conclusions

To conclude, we highlight the need to better understand fundamental processes controlling future EBUS changes under global warming in order to improve our ability to project these changes. This study demonstrates that future EBUS changes are highly complex and each EBUS is governed by a unique set of dynamics, none of which completely fits into a simple dynamical framework, such as the Bakun hypothesis. An important implication from this study is that resolving the core of LLCJs and their rapid drop-off towards upwelling coasts is as important, if not more important, than resolving small-scale ocean features within upwelling zones, because these atmospheric attributes can give rise to the fine-scale near-coast wind stress curls that are vital in simulating and predicting upwelling changes. Thus, robust and reliable projections of future EBUS changes requires high-resolution fully coupled climate models capable of representing both small-scale LLCJ and coastal ocean circulation response to future warming. The regional downscaling approach using eddy-resolving ocean-only models can potentially miss an important part of the dynamics driven by small-scale wind stress curl changes.

Last but not the least, some of the projected changes in EBUSs by HR show competing effects between upwelling intensity and temperature, which may have complicated implications for future changes in nutrient supply and fisheries in EBUSs. For example, surface warming-associated stratification increase can potentially lead to a local reduction in nutrient supply, while an increase in upwelling intensity can enhance nutrient supply. Remote influences from the tropics, possibly driven by variations in the positions of atmospheric pressure systems or large-scale effects caused by e.g., the Pacific Decadal Oscillation or the El Niño-Southern Oscillation42, may produce changes in source water chemistry that affect nutrient supply to the upwelling waters. How these competing changes collectively affect future fisheries environment in EBUSs remains to be investigated43. Future studies need to include ocean biogeochemistry and fisheries components in the high-resolution model framework to assess the full impact of future upwelling changes on marine ecosystems and fisheries.

Methods

CESM simulations

All the simulations38 are based on CESM version 1.3 (CESM1.3) whose atmospheric component is the Community Atmosphere Model version 5 (CAM5) with the Spectral Element dynamical core, the ocean component is the Parallel Ocean Program version 2 (POP2), the sea-ice component is the Community Ice Code version 4 (CICE4), and the land component is the Community Land Model version 4 (CLM4). HR uses a horizontal resolution of ~0.1° for ocean and sea-ice, and ~0.25° for the atmosphere and land, while LR uses a nominal 1° horizontal resolution for all the components. Both HR and LR have a 500-year preindustrial control (PI-CTRL) climate simulation and a historical-and-future transient (HF-TNST) climate simulation from 1850 to 2100, branched from PI-CTRL at year 250, following the protocol for the Coupled Model Intercomparison Project phase 5 (CMIP5) experiments. PI-CTRL was forced by a perpetual climate forcing that corresponds to the 1850 conditions, while HF-TNST was forced by the observed climate forcing until 2005 after which the climate forcing follows RCP8.5 emission scenario. For HR, two additional HF-TNST simulations were performed, branched from the 1850-2100 HF-TNST in 1920 with slightly different atmospheric initial conditions. Combining  the 1850-2100 HF-TNST with the two 1920-2100 HF-TNST simulations gives rise to the three-member ensemble for HR-TNST.

Heat budget analysis

The heat budget analysis for temperature trends in the upwelling region is based on the model temperature equation. First, at each grid point, the temperature equation is vertically integrated from the free surface η to z = 50 m and then integrated in time, which gives

$$\underbrace{\bar{T}(t)-\bar{T}({t}_{0})}_{{{{{{\rm{Temp}}}}}}.\,{{{{{\rm{Change}}}}}}}=\underbrace{-{\int }_{{t}_{0}}^{t}\left(\overline{u\frac{\partial T}{\partial x}}+\overline{v\frac{\partial T}{\partial y}}\right)dt-{\int }_{{t}_{0}}^{t}\overline{w\frac{\partial T}{\partial z}}dt}_{{{{{{\rm{Advection}}}}}}}+\underbrace{{\int }_{{t}_{0}}^{t}\frac{1}{{c}_{P}{\rho }_{0}}({Q}_{net}-SW{|}_{z=-h})dt}_{{{{{{\rm{Net}}}}}}\,{{{{{\rm{Atmos}}}}}}.\,{{{{{\rm{Heating}}}}}}}+\underbrace{R}_{{{{{{\rm{Mixing}}}}}}},$$
(1)

where \(\bar{(\cdot )}\) represents the vertical integral. The term on the left-hand side (LHS) shows changes in temperature over the upper 50 m relative to its initial value at t = t0; the first and second term on the right-hand-side (RHS) are the contribution to the temperature from horizontal and vertical advection of heat, respectively; the third term is the contribution from the net heating in the upper 50 m from the atmosphere, which consists of the net surface heat flux Qnet minus the penetrative solar radiation at z = 50 m; and the last term is the contribution from all turbulence mixing. From (1), we compute the linear trend for each term and then area-average the trends over the 100 km wide strip off each upwelling coast within the green box region shown in Fig. 2a–e. We perform this analysis for ensemble mean of HF-TNST and for PI-CTRL. The results shown in Fig. 3a are the difference between HF-TNST and PI-CTRL, so that any trends due to model drift in PI-CTRL are subtracted out from HF-TNST. We note that in the model the advection term is formulated in a flux form, i.e., ∇⋅vT, in (1) to make the numerical scheme conservative. Here we used a non-flux form, \({{{{{\bf{v}}}}}}\cdot \nabla T\), for the heat budget analysis, because \(\nabla \cdot {{{{{\bf{v}}}}}}T={{{{{\bf{v}}}}}}\cdot \nabla T+T\nabla \cdot {{{{{\bf{v}}}}}}\), and \(\nabla \cdot {{{{{\bf{v}}}}}}\,{{{{{\boldsymbol{=}}}}}}\,0\) due to the incompressibility assumption adopted in POP2. The reason the non-flux form is used is that, even though their sum is zero, each component of \(T\nabla \cdot {{{{{\bf{v}}}}}}\) is non-zero, and these non-zero components do not contribute to the advection process physically. Therefore, one does not want to include them when decomposing the advection into horizontal and vertical components as shown in Fig. 3. We further decomposed the advection into eddy and mean flow induced component, i.e.,

$$\underbrace{\bigg\langle u\frac{\partial T}{\partial x}\bigg\rangle +\bigg\langle v\frac{\partial T}{\partial y}\bigg\rangle +\bigg\langle w\frac{\partial T}{\partial z}\bigg\rangle }_{{{{{{\rm{Total}}}}}}}=\underbrace{\langle u\rangle \frac{\partial \langle T\rangle }{\partial x}+\langle v\rangle \frac{\partial \langle T\rangle }{\partial y}+\langle w\rangle \frac{\partial \langle T\rangle }{\partial z}}_{{{{{{\rm{Mean}}}}}}}+\underbrace{\bigg\langle u^{\prime} \frac{\partial T^{\prime} }{\partial x}\bigg\rangle +\bigg\langle v^{\prime} \frac{\partial T^{\prime} }{\partial y}\bigg\rangle +\bigg\langle w^{\prime} \frac{\partial T^{\prime} }{\partial z}\bigg\rangle }_{{{{{{\rm{Eddy}}}}}}},$$
(2)

where 〈⋅〉 represents monthly mean and \(\left(\cdot \right)^{\prime}\) represents departure from monthly mean. In deriving (2), we again used the incompressibility assumption ∇⋅v  = 0. We output monthly mean 〈uT〉, 〈vT〉, 〈wT〉, 〈u〉, 〈v〉, 〈w〉, and 〈T〉, from which we compute the total advection \(\langle u\frac{\partial T}{\partial x}\rangle +\langle v\frac{\partial T}{\partial y}\rangle +\langle w\frac{\partial T}{\partial z}\rangle =\frac{\partial \langle {uT}\rangle }{\partial x}+\frac{\partial \langle {vT}\rangle }{\partial y}+\frac{\partial \langle {wT}\rangle }{\partial z}\), and the mean advection \(\left\langle u\right\rangle \frac{\partial \left\langle T\right\rangle }{\partial x}+\left\langle v\right\rangle \frac{\partial \left\langle T\right\rangle }{\partial y}+\left\langle w\right\rangle \frac{\partial \left\langle T\right\rangle }{\partial z}\). The eddy advection \(\langle {u{{\hbox{'}}}}\frac{\partial {T{{\hbox{'}}}}}{\partial x}\rangle +\langle {v{{\hbox{'}}}}\frac{\partial {T{{\hbox{'}}}}}{\partial y}\rangle +\langle {w{{\hbox{'}}}}\frac{\partial {T{{\hbox{'}}}}}{\partial z}\rangle\) is then computed as the residual between total and mean advection. We finally estimate the mean horizonal \(\left\langle u\right\rangle \frac{\partial \left\langle T\right\rangle }{\partial x}+\left\langle v\right\rangle \frac{\partial \left\langle T\right\rangle }{\partial y}\,\)and vertical advection \(\left\langle w\right\rangle \frac{\partial \left\langle T\right\rangle }{\partial z}\), respectively.

Vertical velocity at Ekman depth and upwelling index

Ekman depth is computed as the depth at which the simulated cross-shore flow Ucrossshore becomes 20% of the surface Ekman flow UEkman34. The simulated vertical velocity is then interpolated to the Ekman depth at each native model grid, which is then used for computing the upwelling index which is based on the averaged vertical velocity in each coastal bin of 1°(cross-shore) × 1° (alongshore) for LR and ~0.5° (cross-shore) × 0.1° (alongshore) for HR, respectively, to reflect the model resolution difference.