Projected intensification of summer upwelling along the Southeastern Coast of China

Coastal upwelling is a key dynamical process along the southeastern coast of China, regulating regional circulation, biogeochemical cycling, and marine ecosystem dynamics. While previous studies have provided valuable insights into the mechanisms and variability of this coastal upwelling, its future evolution and underlying physical drivers under climate change remain unclear. Here, we investigate projected changes in coastal upwelling along the southeastern coast of China during 2006–2100 using a high-resolution (0.1°) coupled-model realization under the Representative Concentration Pathway 8.5 scenario. During 2006–2100, the area-weighted mean oceanic vertical velocity at 30 m increases at 0.755 m day -1 century -1, while the upwelling area increases at 1.822 × 10³ km² century -1. The trend in the total upwelling volume transport is approximately 0.27 Sv century -1. Enhanced Ekman transport makes the largest contribution (0.26 Sv century -1), accompanied by a geostrophic-transport trend of 0.07 Sv century -1. Marine heatwave (MHW) metrics based on local linear detrended SST are generally lower at the selected coastal sites than at offshore sites, with an 18.3% reduction in MHW days, a 33.2% reduction in mean annual duration, and a 2.4% reduction in mean cumulative intensity. These coastal–offshore contrasts are to some extent consistent with a potential moderating influence of upwelling on MHWs, though rigorous attribution and mechanistic analyses of this effect remain open for dedicated future investigations. This study can advance our understanding of ocean dynamics along the southeastern coast of China under future climate change.

Coastal upwelling prevails along China’s coastline, especially across its southeastern continental margin (Hu and Wang, 2016). During summer, a stable coastal upwelling system develops along the southeastern coast of China under the combined influences of the southwest monsoon, shelf topography, and regional circulation (Gan et al., 2009; Wang et al., 2014). By pumping nutrient-rich subsurface waters into the euphotic zone, coastal upwelling sustains primary marine productivity and shapes regional ocean dynamics and ecosystem functioning (Yin and Harrison, 2008; Yao and Wang, 2021). That said, intensification in upwelling is not unequivocally beneficial. Enhanced organic matter loading can accelerate microbial decomposition and respiratory oxygen depletion, alter coastal oxygen dynamics and potentially exacerbate the risk of localized hypoxia (Chen et al., 2021). Understanding how coastal upwelling along China’s southeastern coast responds to future climate change is therefore essential for assessing the consequences of physical and ecological environmental changes.

Previous studies of coastal upwelling along the southeastern coast of China, especially in the northern South China Sea (NSCS), have largely focused on its formation mechanisms and variability on seasonal to decadal timescales (Zhang et al., 2018; Shu et al., 2018; Zhang, 2026). These studies have demonstrated that the wind stress, shelf topography, background circulation, and large-scale climate variability, including El Niño–Southern Oscillation (ENSO), all contribute to variations in upwelling intensity and spatial extent (Xie et al., 2003; Jing et al., 2011; To Duy et al., 2022). For example, observations and reanalysis products suggest substantial decadal variabilities in the intensity and spatial extent of NSCS upwelling over recent decades, associated in part with changes in regional wind forcing (Liu et al., 2023). In addition to wind-driven processes, nearshore geostrophic adjustment has been identified as an important contributor to coastal upwelling changes (Liao et al., 2024). Despite these advances, most previous studies have focused on historical changes, and the future upwelling evolution along China’s southeastern coast under climate change remains poorly constrained. It remains unclear whether upwelling will intensify or weaken in a warming climate and how changes in wind-driven Ekman transport and geostrophic adjustment will contribute to its future evolution.

Coastal upwelling is not only a crucial dynamical process for associated circulations but also a key mechanism regulating the thermal structure of continental shelf seas (Varela et al., 2015). The upward transport of cold subsurface waters cools the sea surface and limits heat accumulation in coastal regions, potentially modulating the development and persistence of marine heatwaves (MHWs; Santos et al., 2012; Seabra et al., 2019; Varela et al., 2021). However, the global warming is expected to increase the vertical stratification, potentially altering the strength of coastal upwelling and the efficiency of subsurface cold water transport (Cheng et al., 2019), thereby affecting its capacity to regulate MHWs. Although recent studies have demonstrated links between upwelling and the formation and vertical structure of MHWs (Elzahaby and Schaeffer, 2019; Schaeffer et al., 2023; Hu and Wang, 2025), it remains unclear whether future changes of the upwelling will enhance or diminish its influence on coastal MHWs (Wang et al., 2023). This uncertainty is particularly relevant for the southeastern coast of China, where coastal upwelling and extreme ocean warming coexist within a densely populated and economically important marginal sea.

Building upon these gaps, this study uses a high-resolution (0.1°) coupled Earth system model simulation to investigate the projected evolution of coastal upwelling, its dynamic mechanisms, and its underlying association with MHWs along the southeastern coast of China during 2006–2100. The remainder of the paper is organized as follows: Section 2 describes the data and methods; Section 3 presents the main results, including the projected changes of coastal upwelling, major dynamical mechanisms, and its underlying association with MHW characteristics; Section 4 discusses limitations of this study and implications of the results; and Section 5 summarizes the main conclusions.

Model simulations

To investigate climate change impacts on the upwelling system along the southeastern coast of China, we employed an eddy-resolving global climate simulation from the International Laboratory for High-Resolution Earth System Prediction (iHESP; Chang et al., 2020; Zhang et al., 2020). Built upon the Community Earth System Model (CESM) framework, the iHESP coupled model has been widely applied in prior coastal upwelling studies (Jing et al., 2023; Liao et al., 2024). Its ocean component employs a curvilinear grid (Supplementary Figure S1) with a nominal 0.1° horizontal resolution and 62 vertical layers (Chang et al., 2020; Zhang et al., 2020). Compared with most coarse-resolution (0.25~1°) global climate models, the 0.1° iHESP configuration provides a relatively detailed representation of coastal geometry and shelf bathymetry (Hewitt et al., 2020; Jing et al., 2023) and can to some extent provide useful insights into regional upwelling research (Jing et al., 2009). However, nominal grid spacing of 0.1° does not guarantee that narrow coastal upwelling fronts and the shallow thermocline are fully resolved. Supplementary Figure S1 shows bathymetry and native model grid. Overall, the bathymetry and four cross-shelf sections analyzed can be reasonably represented in iHESP, although finer processes (e.g., submesoscale) are largely misrepresented.

We used model outputs under the Representative Concentration Pathway 8.5 (RCP8.5) scenario from 2006 to 2100. Monthly mean ocean temperature, salinity, and velocity were analyzed to investigate changes in upwelling intensity, spatial extent, volume transport and the associated thermal state changes. Barotropic geostrophic currents are derived from sea surface height data, while baroclinic adjustment processes are examined via thermohaline structure analysis. Wind stress and wind-stress curl are used to compute Ekman transport and pumping, allowing us to evaluate the contribution of wind-driven dynamics. Additionally, daily sea surface temperature (SST) is used to identify MHW events, facilitating the analysis of the spatial association between the coastal upwelling and MHW characteristics. Our temporal analysis compares the historical reference period (P1: 2006–2025) and the late-century period (P2: 2081–2100).

To further evaluate the model’s regional applicability, we compared simulated summer SST with satellite-based SST from the Operational Sea Surface Temperature and Ice Analysis (OSTIA) over 2006–2025 (Supplementary Figure S2). The model reproduces the broad coastal–offshore temperature contrast, including relatively cool waters along the Guangdong–Fujian coast and the western Taiwan Strait. Regional mean summer SST is 27.17 °C in iHESP and 27.49 °C in OSTIA, with a regional root mean square error (RMSE) of 0.81 °C. Given that SST serves as a well-established proxy for coastal upwelling along China’s coast, this comparison bolsters our confidence in applying iHESP to upwelling analyses.

The study region and upwelling diagnostics

Our study region covers the southeastern coast of China, mainly encompassing the Zhejiang–Fujian–Guangdong coastal waters. As shown in Figure 1A, the southwest monsoon prevails over the study region annually from June to August (i.e., boreal summer). Based on the simulated climatological summer mean vertical velocity, coastal upwelling is primarily concentrated off eastern Guangdong and Fujian (Figure 1B), with weaker upwelling along the Zhejiang coast. We analyze the oceanic vertical velocity at 30 m to characterize summer coastal upwelling. This depth is selected based on previous observations showing distinct upwelling-related temperature and salinity signatures at 30 m along the Fujian–Guangdong coast (Pan et al., 2012; Wang et al., 2016).

Within the fixed study polygon, we identify summer upwelling as wet grid cells (magenta polygon in Figure 1B) where the summer-mean vertical velocity at 30m exceeds 1.0 m day -1. The threshold of 1.0 m day -1 is adopted to filter noisy patterns. A sensitivity analysis of this threshold was performed, and the main results and conclusions remained robust when the threshold varied between 0.5 and 1.5 m day -1 (Supplementary Figure S3). Upwelling intensity is the area-weighted mean vertical velocity over these cells, and upwelling area is the sum of their areas.

Volume flux budget diagnosis

To quantitatively analyze the dynamical mechanisms underlying future changes in upwelling along the southeastern coast of China, we calculated volume fluxes following Liao et al. (2024). For a coastal region with open lateral boundaries (Figure 2), the principle of volume conservation (the seawater is recognized as incompressible) dictates that the vertical upward volume flux through a bottom reference layer equals the net horizontal outflow flux across the lateral open boundaries. We further decomposed the lateral horizontal divergence into wind-driven Ekman transport and geostrophic transport (Jing et al., 2023), as expressed in Equation 1:

where V U p represents the total vertical volume transport through the horizontal reference surface at the selected depth (Figure 2), with positive values indicating upward motion. V E k represents the net outward Ekman transport across the lateral boundaries, diagnosed from the area-integrated horizontal divergence of Ekman transport. This diagnostic accounts for both coastal divergences associated with cross-shelf transport associated with wind stress and Ekman pumping associated with wind-stress curl (Jing et al., 2023). V G e o represents net outflow across the open boundaries due to lateral geostrophic currents, with contributions associated with sea-surface-height gradients and horizontal density variations. Positive values of these lateral transport components indicate net outflow, whereas negative values indicate net inflow.

For the selected polygonal area, total vertical transport is calculated by integrating the vertical velocity over the selected horizontal area (Equation 2):

where S is the upwelling area, and W is vertical velocity. dx, dy are the grid spacings in the zonal and meridional directions.

Following Jing et al. (2023), V E k is calculated by integrating the horizontal divergence of Ekman transport over the study region (Equation 3):

where U E and V E are the horizontal Ekman transports in the zonal and meridional directions, respectively, and are calculated using Equations 4 and 5:

where τ x and τ y represent the zonal and meridional wind stresses, respectively. ρ 0 is the reference density of seawater, and f is the Coriolis parameter. V G e o can be calculated via vertical and horizontal integration along the polygon’s open boundaries based on the horizontal geostrophic velocity u g, v g of each layer (Equation 6; Jing et al., 2023; Liao et al., 2024):

where u g and v g are the geostrophic velocity components in the zonal and meridional directions, respectively.

We also calculated the squared buoyancy frequency (N 2) to assess changes in upper-ocean density stratification (Equation 7):

where g is gravitational acceleration and ρ is potential density. Vertical potential-density gradients were computed using central differences on the model’s vertical grid. Larger values of N 2 indicate stronger density stratification, which tends to inhibit vertical mixing (Li et al., 2020).

Detection of MHWs

MHWs are defined as discrete extreme warming events where SST exceeds the local climatological 90th percentile threshold for at least five consecutive days. Successive events separated by two days or less are combined into a single MHW (Hobday et al., 2016). Traditional MHW detection relies on fixed historical baselines in which the climatological mean and 90th-percentile threshold are calculated from the historical reference period and then held constant when identifying subsequent events (Hobday et al., 2016). For a baseline comparison (Supplementary Figure S4), we applied the conventional fixed-baseline method to undetrended SST, using 2006–2025 to calculate the climatological mean and threshold (Equation 8):

where S S T (t) denotes daily SST at time t. For calendar day j, I j contains the sampling dates within an 11-day window centered on that day in each year of the reference period [ t s t a r t, t e n d ] and N j is the number of samples. S S T m (j) and S S T 90 (j) are the climatological mean and 90th-percentile threshold, respectively. P 90 denotes the 90th percentile.

However, in century-scale climate projections, persistent background warming could lead to a future state where marine heatwave conditions prevail across nearly the whole year (Amaya et al., 2023; Rosselló et al., 2023), a condition referred to as a “perpetual heatwave state” (e.g., Supplementary Figure S4). Such prolonged exceedance makes it difficult to distinguish changes in temperature variability from the effects of long-term mean warming (Frölicher et al., 2018; Wang et al., 2023).

Consequently, we implemented a partial moving baseline method incorporating local linear detrending (LLD, Wang et al., 2022). This method involves two separate detrending steps. First, a single linear fit was estimated from daily SST over the reference period (2006–2025) and subtracted from that series. The detrended reference SST series was then used to calculate the climatological mean and threshold following the procedure above. Then, for each target year, a linear fit was estimated from daily SST within an associated 31-year window and subtracted from the target-year SST. The window was centered on the target year where possible and shifted near the record boundaries to retain its 31-year length. The detrended target-year series were concatenated, and MHWs were identified using the detrended target-year SST (Equation 9) and the threshold calculated from the detrended reference SST series (Equation 10).

where S S T L L D (t) represents the target year SST after detrending. β w and α w are the slope and intercept of the linear fit estimated over the 31-year window associated with the target year. S S T r e f (t) is the reference-period SST after subtraction of a single linear fit estimated exclusively over 2006–2025, with slope β f and intercept α f.

For spatial sampling, we selected valid grid cells satisfying two simultaneous criteria: SST below the median (S S T < S S T 50 t h) and positive vertical velocity (W > 0). Among these cells, we selected six grid points with the highest interannual variability in vertical velocity as the coastal upwelling sites (Figure 1C). For each core site, we selected one offshore control point at the same latitude, extending approximately 40 longitudinal grids (about 4°) eastward, where summer upwelling dynamics are climatologically absent. These six coastal–offshore pairs are used to describe spatial contrasts in total MHW days, duration, and cumulative intensity. When calculating multi-year summary statistics, years without MHW events are included, with their annual MHW metrics assigned a value of zero.

Supplementary Figure S4 shows that annual total MHW days approach year-round occurrence in the late-century period when adopting the fixed historical baseline, whereas MHW days derived from LLD remain well below this level. We therefore adopted the LLD method for subsequent MHW analyses to examine extremes after removing the estimated local warming trend.

Decomposition of vertical temperature advection

To examine how changes in upwelling and thermal structure jointly affect the vertical advective cooling tendency, we calculated a vertical temperature-advection diagnostic, denoted by Q u p w. It was defined at the selected depth as:

where T is temperature, and ∂ T / ∂ z represents vertical temperature gradient. Here, vertical velocity and the vertical temperature gradient determine the local temperature tendency due to vertical advection. Negative Q u p w indicates a cooling tendency. It should be noted that Equation 11 does not represent the full upper-ocean heat budget, which also includes surface heat fluxes, horizontal advection, entrainment, mixing, and diffusion.

To separate the contributions of changes in vertical velocity and temperature structure, we expressed the late-century fields (P2: 2081–2100) as the sum of the reference-period means (P1: 2006–2025) and perturbations induced by future climate change (i.e., P2 minus P1 changes):

Substituting Equations 12 and 13 into Equation 11, the change in vertical temperature advection can be decomposed into three algebraic terms:

where − Δ W · (∂ T / ∂ z ¯) P 1 is the contribution from changes in vertical velocity while holding the P1 temperature gradient fixed, − W ¯ P 1 · Δ (∂ T / ∂ z) is the contribution from changes in the vertical temperature gradient while holding the P1 vertical velocity fixed, and − Δ W · Δ (∂ T / ∂ z) represents the nonlinear term arising from simultaneous changes in both quantities. These terms describe the algebraic contributions to the simulated change in local vertical temperature advection, rather than independent causal contributions to the complete upper-ocean heat budget.

Projected changes in coastal upwelling

During 2006–2100, summer upwelling along the southeastern coast of China exhibits increasing trends in both intensity and threshold-defined area in the analyzed RCP8.5 simulation (Figure 3). The area-weighted mean vertical velocity at 30 m increases at a rate of 0.755 m day -1 century -1. Mean summer upwelling intensity in P2 is 7.8% higher than that in P1. Concurrently, the upwelling area increases at a rate of 1.822×10 3 km² century -1. These results demonstrate the intensification and expansion of the threshold-defined upwelling region within this simulation.

Concordant with the enhanced upwelling, the coastal summer wind stress and wind-stress curl also exhibit significant positive trends (Figures 3C, D), increasing at rates of 0.0145 N m - ² century -1 and 1.995×10 -7 N m - ³ century -1, respectively. Mean summer wind stress and wind-stress curl in P2 are 45.5% and 54.3% higher than in P1, respectively. The synchronous enhancement of the wind field indicates a strengthening of the summer monsoon forcing, which provides a more favorable dynamical background for coastal Ekman transport and pumping. We also examine the accompanying spatial pattern of ocean warming, which is heterogeneous across this marginal sea (Supplementary Figure S5). Warming rates in the surface layer and outer shelf generally exceed 3.5 °C century -1, whereas warming is weaker in parts of the middle and lower water column along the upwelling pathways. In particular, regions of weaker warming at approximately 25–60 m along T1 and T4 spatially coincide with the simulated upwelling cores (Supplementary Figures S5, S6).

To further examine the vertical structure of the upwelling changes, we compared summer vertical velocities along four cross-shelf transects between P1 and P2 (Figure 4). In these transects, upward motion is concentrated mainly at depths of approximately 20–50 m, with vertical-velocity maxima generally occurring near 30 m, therefore supporting the selection of 30 m as the reference depth for the upwelling diagnostics (Figures 1B, 3). Compared to P1, the maximum vertical velocities along transects T1–T4 increase by 2.9%, 8.5%, 6.2% and 17% in P2, respectively. The enhancement occurs mainly in shallow coastal waters and along the sloping shelves. Thus, the transect maxima consistently increase, although the magnitudes and spatial patterns of change differ among transects.

Mechanism analysis of the upwelling intensification

To evaluate the dynamical drivers of future upwelling intensification along the southeastern coast of China, we decompose the regional total vertical transport into wind-driven Ekman and geostrophic components (Figure 5). Results indicate that the long-term trend in the total regional vertical transport is approximately 0.27 Sv century -1. The largest diagnosed component is the wind-driven Ekman transport trend, approximately 0.26 Sv century -1. Geostrophic transport adjustments provide a secondary and much weaker contribution, increasing by about 0.07 Sv century -1.

The Ekman and geostrophic transport trends sum to approximately 0.33 Sv century -1, exceeding the diagnosed total vertical transport trend of approximately 0.27 Sv century -1. The residual, defined as total vertical transport minus the Ekman and geostrophic contributions, has a trend of approximately −0.06 Sv century -1. This residual may arise from several sources, including unresolved ageostrophic non-Ekman lateral transport, differences in spatial discretization among the diagnostics, and inconsistencies in depth integration near shallow coastal boundaries. Monthly output may also limit the precise diagnosis of these processes.

Time evolution further shows that Ekman transport and total vertical transport broadly co-vary (Figure 5). Between P1 and P2, their mean values increase by 0.16 and 0.17 Sv, respectively, while mean geostrophic transport increases by 0.07 Sv. These results identify wind-driven Ekman transport as the largest contribution to the upwelling increase, alongside a smaller positive contribution from geostrophic adjustment.

The vertical profiles show that the increase in total vertical transport (i.e., upwelling intensification) is concentrated mainly at depths of 20–50 m, with negligible change below 100 m (Figure 5B). However, changes in the geostrophic contribution are depth dependent. Geostrophic transport above 50 m increases by approximately 0.02–0.1 Sv in P2 relative to P1, whereas transport at 50–200 m decreases by approximately 0.05–0.14 Sv. Thus, the simulated upwelling intensification does not reflect a uniform strengthening of dynamics across the entire water column, but rather a wind-driven surface intensification accompanied by a reorganization of the regional ocean circulation.

We further partitioned the total geostrophic transport into barotropic and baroclinic components (Supplementary Figures S7, S8). The diagnosed increase is dominated by the barotropic component (94%), with a much smaller baroclinic contribution (6%). The nearshore–offshore SSH contrast increases between P1 and P2 (Supplementary Figure S7). This enhanced cross-shelf SSH gradient is consistent with the geostrophic adjustment. Cross-shelf barotropic gradients along T1–T4 increase by 6.2%–13.2%, with the largest increase along T4. The N 2 diagnostics show local increases and decreases between P1 and P2, rather than a uniform strengthening of upper-ocean density stratification (Supplementary Figure S8).

Coastal upwelling and MHWs under future warming

Coastal upwelling may moderate marine heatwaves by bringing cooler subsurface water into the upper ocean (Izquierdo et al., 2022). To examine these potential effects, we compared LLD-based MHW characteristics at the selected coastal upwelling sites and their offshore counterparts (Figure 1C), with the results shown in Figure 6 and Table 1. We further analyzed changes in local vertical temperature advection between P1 and P2 to examine how evolving upwelling and thermal structure jointly affect the local advective cooling tendency (Figure 7), providing physical context for the potential influence of upwelling on MHW changes.

Mean annual total MHW days at the coastal sites are 5.97 days fewer than at the offshore sites (Figure 6B; Table 1). With no-event years included as zeros, the mean annual duration is 9.34 days at the coastal sites and 13.98 days offshore, a difference of 33.2% (Figure 6D; Table 1), and the mean cumulative intensity is 2.4% weaker at the upwelling sites (Figure 6F; Table 1). These coastal–offshore contrasts are consistent with a potential moderating influence of upwelling on MHWs.

Notably, under the LLD-based threshold definition, total MHW days exhibit a downward trend in both regions, with offshore waters declining faster (-1.20 days decade -1) than the coastal upwelling sites (-0.54 days decade -1; Figure 6A). These trends describe MHW changes after removal of the estimated local warming trend and do not imply declining absolute SST or reduced absolute thermal exposure. Therefore, the long-term trend in MHW days alone cannot represent the change in true MHW risk. In contrast, duration and cumulative intensity provide complementary measures of the persistence and accumulated thermal anomaly of MHW events. Projections indicate that MHW duration in the open ocean will increase (0.15 days decade -1), while that in the upwelling sites will decrease (-0.25 days decade -1; Figure 6C). Concurrently, the declining cumulative intensity in the upwelling sites (-0.55 °C·day decade -1) is associated with reduced accumulated thermal anomalies (Figure 6E), although it does not imply the elimination of MHW occurrence.

To examine the accompanying changes in the local thermal effect of upwelling, we decomposed the P2-minus-P1 changes of the vertical temperature advection into three components (Figure 7; Equation 14). The magnitude of the local vertical temperature advection increases by approximately 0.026 °C day -1 (Figures 7A, B). The vertical-velocity and vertical-temperature-gradient terms account for approximately 53.3% and 41.9% of this diagnosed change, respectively, while their interaction accounts for the remaining 4.8%. The velocity term reflects the contribution from changing vertical velocity while retaining the P1 temperature gradient. The gradient term reflects the contribution from changing vertical thermal structure while retaining the P1 vertical velocity.

These results show that the simulated increase in local advective cooling results from changes in both upward transport and the vertical temperature gradient (Figure 7). Enhanced upwelling improves the capacity to transport subsurface cold water to the surface, and the strengthened vertical temperature gradient further increases the local thermal contribution of this vertical transport. The interaction term represents the product of their concurrent changes (Figure 7B), rather than independent evidence of physical feedback. This decomposition captures the local vertical-advection tendency and merely demonstrates the potential of upwelling to moderate MHWs. In other words, it does not quantify the full upper-ocean cooling response or directly attribute the coastal–offshore MHW differences to upwelling.

Coastal upwelling in monsoon-dominated marginal seas is closely linked to local wind forcing through Ekman transport (Xie et al., 2003; Wang et al., 2014; Hu and Wang, 2016), while the concurrent geostrophic changes indicate accompanying adjustments in regional circulation (Figure 5). The residual in the transport budget nevertheless limits precise attribution of the projected change to individual components. Upwelling also affects local temperature through vertical temperature advection, whose cooling tendency depends on both vertical velocity and the vertical temperature gradient (Figure 7). Considering these quantities together therefore helps explain how changes in circulation and thermal structure jointly shape the thermal effect of upwelling, providing physical context for its potential influence on MHWs.

The interpretation of the results in this study is subject to uncertainties associated with model resolution and experimental coverage. The nominal 0.1° ocean grid is finer than the 1/6° grid used in an earlier regional upwelling study (Jing et al., 2009), but narrow coastal fronts, sharp vertical temperature gradients, and localized vertical-velocity maxima may remain insufficiently resolved. These limitations may affect the magnitude of the temperature-advection diagnostic and its decomposed contributions. The comparison with OSTIA supports the model’s representation of broad coastal thermal patterns, although it does not directly validate the subsurface dynamics. Furthermore, the analysis uses a single realization under RCP8.5 and therefore cannot independently distinguish the externally forced response from internal variability or quantify uncertainty across models and scenarios. The coupled simulation provides a useful framework for examining regional changes, while additional realizations, models, scenarios, and finer-resolution regional experiments are needed to assess their robustness.

Interpretation of the potential upwelling–MHW relationship also requires distinguishing the reference conditions used to identify MHWs. The fixed historical baseline retains the influence of background warming, leading to prolonged future threshold exceedance, whereas the LLD method identifies extremes after removal of the estimated local warming trend (Supplementary Figure S4). Under LLD, annual MHW days decline at both coastal and offshore sites, more rapidly offshore (Figure 6A), but this does not imply reduced thermal exposure relative to historical conditions.

The selected coastal upwelling sites generally exhibit lower LLD-based MHW statistics than their offshore counterparts (Figure 6; Table 1). These contrasts are consistent with a potential moderating influence of upwelling, but the sampling design does not isolate that influence. The coastal sites were selected for their high interannual variability in vertical velocity, and coastal–offshore differences in bathymetry, thermal structure, circulation, and atmospheric forcing may also contribute to the contrasts. Accordingly, the spatial comparison does not establish that future upwelling intensification causes weaker MHWs. The vertical temperature-advection diagnostic provides complementary physical information, but a complete assessment of MHW temperature changes would also require surface heat fluxes, horizontal advection, and mixing.

Furthermore, it is worthwhile to note that upwelling and thermal structure may interact in a two-way feedback loop. Vertical transport of upwelling modifies temperature and density, while the resulting structure can influence circulation, mixing, and the thermal effect of subsequent transport. Such interactions may influence the persistence and thermal expression of upwelling and its potential modulation of MHWs. However, the interaction term in our decomposition represents the product of concurrent changes in vertical velocity and the temperature gradient; it does not independently demonstrate a feedback loop. Targeted process experiments and time-resolved heat-budget analyses would help establish the direction and relative importance of these interactions and clarify their consequences for marine heat extremes.

Using the single iHESP coupled-model realization with a nominal ocean resolution of 0.1° under the RCP8.5 scenario, we investigated the projected changes in coastal upwelling along the southeastern coast of China during 2006–2100, their dynamical mechanisms, and potential links with MHW changes.

Within the analyzed simulation, the area-weighted mean oceanic vertical velocity at 30 m increases at a rate of 0.755 m day -1 century -1. The upwelling area increases at a rate of 1.822 × 10³ km² century -1. These results indicate long-term intensification and expansion of the coastal upwelling region along China’s southeastern coast under the analyzed high-emission scenario.

Wind-driven Ekman transport is the largest diagnosed component of the vertical-transport trend, increasing at approximately 0.26 Sv century -1, while geostrophic transport increases at approximately 0.07 Sv century -1. The concurrent changes indicate a wind-driven response accompanied by regional circulation adjustment. However, the sum of Ekman and geostrophic contributions exceeds the total vertical-transport trend, with a residual trend of approximately −0.06 Sv century -1. The incomplete budget closure limits precise attribution of the total change.

The baseline comparison highlights the distinction between MHWs defined relative to historical thermal conditions and those identified after removal of the estimated local warming trend. Using the LLD method and the averaging convention, the selected coastal upwelling sites exhibit 18.3% fewer annual MHW days, a 33.2% lower mean annual duration statistic, and 2.4% lower mean annual cumulative intensity than the offshore comparison sites. The vertical temperature-advection diagnostic further indicates that changes in both vertical velocity and the vertical temperature gradient contribute to a stronger cooling tendency. The results are consistent with a potential moderating influence of upwelling on MHW characteristics, although the present analysis does not establish causality.

This study provides a regional assessment of coastal upwelling changes along the southeastern coast of China under future warming. It also provides a physical basis for further investigations on how upwelling environments may modulate marine heat extremes under future warming. Combining these diagnostics with ecosystem and biogeochemical modeling could help assess the consequences of changing upwelling for coastal productivity and ecosystem vulnerability in the future.