Sea Ice Loss leads to regime shifts in the arctic biological pump

Domain: earth_science
Published: 2025-11-24
Editorial outcome: minor_revision

Abstract
The Arctic Ocean has undergone accelerated warming and a marked decline in sea ice over recent decades. Yet, the response of the biological pump—a critical mechanism for atmospheric carbon sequestration—remains poorly understood. Here, we develop a satellite-derived dataset (2003–2022) to identify a regime shift in the Arctic biological pump. Between 2003 and 2012, the strength and efficiency of the biological pump increased rapidly, primarily driven by sea ice decline. However, from 2013 to 2022, this trend plateaued, coinciding with stabilized sea ice conditions and increased phytoplankton biomass. Earth system model simulations (1850–2100) support the observed link between biological pump enhancement and sea ice loss, and project a future regime shift towards a weakened biological pump under nearly ice-free conditions, associated with shifts in phytoplankton community structure. These findings underscore the Arctic’s vulnerability to climate-driven changes, with far-reaching implications for Arctic carbon sequestration and ecosystem stability.

Introduction
Arctic sea ice decline represents a striking indicator of global climate change, driving significant alterations in ecosystem processes, including the production and deep-ocean export of organic carbon, or, as it is popularly known, carbon sequestration. The Arctic Ocean (AO) serves as a significant sink for atmospheric carbon dioxide (CO2), accounting for up to 12% of global oceanic carbon uptake, despite comprising less than 5% of the global sea surface area. Much of this CO2 uptake is driven by the biological pump, which transfers carbon from the surface ocean and its linked atmosphere to the deep waters. While the strength and efficiency of the biological pump strongly influence the degree of carbon export to the deep ocean, its response to rapid sea ice loss in the AO remains poorly understood. In particular, Arctic sea ice extent has drastically declined in recent decades, characterized by earlier seasonal retreats and reduced ice thickness. However, the impacts of sea ice decline on the biological pump remain debated. Some studies suggest that sea ice decline reduces vertical carbon export in the AO due to increased stratification and diminished nutrient supply, while others argue that sea ice retreat enhances the biological pump’s strength through improved light penetration, an extended growing season, and intensified ocean circulation. A key reason for these conflicting findings is that observational datasets predominantly consist of temporally sparse measurements, offering only limited coverage of the dynamic ocean environment.

Resolving these conflicting statements has proven challenging due to the region’s spatial and temporal variability and the scarcity of consistent observations. Ocean color satellite-based observations have been utilized to continuously monitor the variations in the biological pump by evaluating net primary production (NPP). However, this method depends on empirical links between the export efficiency (e-ratio) and satellite-derived NPP and sea surface temperature (SST), which vary regionally across the global oceans, leading to discrepancies in carbon export estimates. Alternatively, Earth system models (ESMs) have been employed to monitor alterations in the biological pump under climate change conditions; however, in most ESMs, carbon export is parameterized almost exclusively via particle sinking, whereas other mechanisms, such as zooplankton vertical migration and the subduction or vertical redistribution of dissolved and suspended organic carbon, are largely omitted. Therefore, a fundamental question remains unanswered: has the biological pump in the AO strengthened or weakened in recent decades as a result of diminishing sea ice extent?

To address this critical knowledge gap in AO carbon cycle research, we present a satellite-based dataset of net community production (NCP) spanning 2003–2022. NCP, representing the net accumulation of organic carbon in the surface ocean after accounting for heterotrophic respiration by consumers, effectively reflects the export flux of carbon from the surface under steady-state conditions and across large spatial scales. A portion of the exported carbon can be transferred to the deep ocean and stored there, mediated by several biological pump processes such as gravitational settling of particles, zooplankton-driven vertical transport, and the subduction or mixing of organic material. Accordingly, NCP serves as a key proxy for biological pump strength. Normalized by NPP, the e-ratio (NCP/NPP) provides an index of carbon export efficiency. Although ship-based measurements have offered foundational insights, their limited spatiotemporal coverage has hindered robust assessments of long-term trends, particularly in the AO, where no decadal-scale, high-resolution monitoring existed prior to this study. Here, we overcome these limitations using a satellite-based framework, enabling systematic quantification of spatial variability and multi-decadal changes in the Arctic biological pump.

In this work, we estimate NCP using a semi-analytical algorithm that integrates satellite observations, reanalysis products, machine learning techniques, and biogeochemical models (see Methods; Supplementary Figs. 1–7, Supplementary Tables 1–2). Additionally, we use NPP to assess variations in the e-ratio. We employ five established models to map NPP distribution, including VGPM, Eppley-VGPM, CbPM, CAFE, and a locally calibrated Arctic Ocean Productivity Model (AOPM) (Supplementary Figs. 8, 9). Our analysis focuses on the period from late spring to late summer (June–September) because sea ice typically begins to melt in early June and refreezes in October for the majority of the AO, and the mixed layer remains relatively stable during this period. We statistically analyze trends in NCP, NPP, and the e-ratio over the past two decades across the entire AO and its 11 sub-regions (Fig. 1a) in relation to sea ice dynamics. We incorporate ESM simulations to extend the analysis from 1850 to 2100, capturing the evolution of the impacts of warming and sea ice loss on the biological pump. Collectively, these datasets provide critical insights into the historical trajectory and the future response of the Arctic biological pump under ongoing climate change.

Results
On average, the AO exhibits a mean NCP of 8.81 ± 2.03 mmol C m−2 d−1, surpassing the upper limit of NCP observed in other open ocean regions (Supplementary Table 3) and highlighting its significant contribution to carbon cycling in the global ocean. The area-integrated NCP is approximately 55.24 ± 15.27 Tg C yr−1, consistent with a previous estimate of 61 ± 23 Tg C yr−1 within the mixed layer depth (MLD). The mean e-ratio across the AO is 0.12, indicating that about 12% of the NPP is exported out of the MLD through the biological pump.

Both NCP and the e-ratio exhibit pronounced subregional heterogeneity (Fig. 1b, c and Supplementary Fig. 10). The highest NCP values (9.62 ± 2.12 mmol C m−2 d−1) occur in inflow shelves, driven mainly by enhanced phytoplankton biomass stimulated by nutrient-rich subarctic inflows. Outflow (7.83 ± 2.83 mmol C m−2 d−1) and interior shelves (6.10 ± 2.69 mmol C m−2 d−1) also exhibit relatively high NCP. Conversely, notably low NCP values appear in river-dominated coastal zones on interior shelves (e.g., the Kara and Laptev Seas), primarily due to substantial inputs of terrigenous dissolved organic matter, most of which is remineralized to dissolved inorganic carbon (DIC) in coastal areas. The e-ratio shows a broadly similar spatial pattern (Fig. 1c). The highest e-ratio values, exceeding 0.5, occur in inflow regions such as the Northeast Barents Sea and Chukchi Sea shelves, whereas the lowest e-ratio values (<0.05 or 5%; Fig. 1e) are observed in interior shelf coastal zones. These low values reflect high NPP combined with intense respiration, fueled by terrestrial organic carbon inputs and rapid organic matter turnover in the shallow shelf (Fig. 1b, c and Supplementary Fig. 9). Interestingly, the central basin, despite its low NPP (Supplementary Fig. 9b), shows moderate NCP (5.26 ± 4.40 mmol C m−2 d−1) and relatively high e-ratio values (reaching up to 0.40, with a mean of 0.13). These patterns are likely driven by phytoplankton blooms near the ice edge, highlighting the pivotal role of the marginal ice zone in deep-ocean carbon sequestration.

Satellite observations over the past two decades reveal distinct trends in the strength of the Arctic biological pump (Fig. 2a, d), with rapid intensification during the first decade (2003–2012) followed by a marked slowdown in the second decade (2013–2022). During the first decade, NCP increased at 0.067 mol C m−2 yr−2, three times the second decade’s rate of 0.021 mol C m−2 yr−2. This acceleration was strongest in the Kara and Barents Seas, where NCP rose by 0.093 and 0.076 mol C m−2 yr−2, respectively (Fig. 2b and Supplementary Fig. 10). Moderate but significant increases also occurred in the Norwegian, East Siberian, and Greenland Seas. In contrast, other subregions, including the Laptev and Beaufort Seas, showed no significant trend (P > 0.05). In 2013–2022, only the Kara Sea maintained a significant increase, while NCP growth elsewhere stalled or reversed (Fig. 2c and Supplementary Fig. 10).

The e-ratio generally tracked NCP variations across the AO, independent of the choice of NPP products (Supplementary Fig. 11). It increased widely in the first decade (0.0117 yr−1 overall) but stabilized in the second decade (Fig. 2d). Early increases were concentrated in the Kara (0.0147 yr−1), Barents (0.0135 yr−1), East Siberian (0.0068 yr−1), and Greenland (0.0132 yr−1) Seas, and detectable—though often not significant—in most other shelves (Fig. 2e and Supplementary Fig. 10). By contrast, in 2013–2022 localized gains occurred only in the Kara Sea and parts of the central basin’s marginal ice zone, while previously increasing regions such as the East Siberian, Chukchi, and Beaufort Seas declined (Fig. 2f and Supplementary Fig. 10).

By controlling the availability of light, nutrients, and habitat space essential for primary producers in the AO, sea ice plays a pivotal role in ecosystem functioning. Accordingly, changes in sea ice exert profound influences that may restructure Arctic marine ecosystems. Between 2003 and 2022, Sea Ice Concentration (SIC) declined significantly while open water duration (OWD) increased (Fig. 3a and Supplementary Fig. 12a), coinciding with NCP changes. From 2003 to 2012, SIC declined rapidly, and NCP increased significantly, particularly in interior and outflow shelves such as the Kara, Laptev, and northern Barents Seas, as well as the Beaufort Sea (Fig. 3b and Supplementary Fig. 12b). In contrast, from 2013 to 2022, the rate of SIC loss slowed, leaving concentrations relatively stable. In this period, SIC decreased slightly—but not significantly—in the eastern interior shelves (East Siberian, Laptev, and Kara Seas), while a slight increase was observed in the western shelves (Beaufort Sea) (Fig. 3c and Supplementary Fig. 12c). Correspondingly, the biological pump remained stable (Fig. 3d).

Correlation analysis furthers supports the coupling between biological pump variations and sea ice changes over the past two decades. A robust negative correlation was identified between NCP and SIC variation (r = −0.72, P < 0.001) during 2003–2022. This correlation was predominantly driven by the strong association between NCP and SIC during the first decade. For instance, NCP was significantly correlated with SIC (r = −0.77, P < 0.01) and OWD (r = 0.81, P < 0.01) in the first decade (Fig. 3d and Supplementary Fig. 12d); however, these correlations were not significant in the second decade. In contrast, during 2013–2022, a significant positive correlation emerged between NCP and chlorophyll-a (Chla) (r = 0.77, P < 0.01, Fig. 3e). The trend in the e-ratio closely followed NCP variation, exhibiting significant correlations with SIC (r = −0.77, P < 0.01) and OWD (r = 0.81, P < 0.01) during 2003–2012. However, no significant correlations were observed between the e-ratio and either sea ice variation or Chla during the second decade. These findings indicate a shift in the primary drivers of the biological pump over the past two decades.

Mechanistic evidence supports the observed coupling between declining SIC and increasing NCP during the early melting period. The expansion of open water areas caused by reduced SIC, combined with more frequent high-latitude storms, markedly weakened sea surface stratification. This transition shifted the upper ocean from a meltwater-stratified regime—driven by surface accumulation of low-salinity meltwater, strong stratification, a shallow MLD, and limited vertical exchange under higher SIC—to a mixed-layer-homogenized regime under lower SIC, where wind-driven turbulence disrupts stratification, deepens the MLD, and enhances vertical mixing. Reduced SIC is also linked to intensified eddy activity along shelf break regions, which facilitates the downward transport of organic matter. Together, these physical changes create favorable conditions for organic carbon export to the deep ocean.

Significant correlations between SIC or OWD and MLD (Fig. 3f; Supplementary Figs. 12e and 13d), along with a positive correlation between NCP or e-ratio and MLD during the active melting period of 2003–2012 (Fig. 3g; Supplementary Figs. 13e and 14), support the hypothesis that reduced upper ocean stability, reflected by MLD deepening under declining SIC, was the key physical driver of enhanced carbon export (higher NCP and higher e-ratio). In the shallow mixed layer of the AO (typically 5–30 m in summer), even small MLD changes can substantially strengthen the biological pump and carbon export efficiency through multiple pathways. For example, MLD deepening entrains nutrient-rich bottom waters, alleviating nutrient limitations in strongly stratified surface layers (see Supplementary Fig. 15). Nutrient-deficient phytoplankton exhibit disproportionate productivity gains from small nutrient inputs. Reduced SIC also markedly increases light penetration into the water column. Deeper mixing redistributes phytoplankton out of the surface light-inhibition zone and optimizes photosynthetic efficiency under enhanced light conditions. Enhanced nutrient supply, combined with optimized photosynthesis, supports higher primary production in marginal shelves and drives a shift in phytoplankton communities toward large-celled taxa such as diatoms. Larger phytoplankton enhance carbon export because their faster sinking accelerates POC transfer to depth. Rapid sinking shortens particle residence time in the euphotic zone and limits remineralization by heterotrophs. As a result, a larger fraction of primary production was transferred to the deep ocean as sinking particles (i.e., the e-ratio; Fig. 2d), contributing to higher NCP. The link between phytoplankton size structure and carbon export explains the rise in NCP, despite no significant increase in NPP during this period (Supplementary Fig. 9d). Together, these findings provide mechanistic, spatial, and statistical evidence that sea ice decline was the primary driver of the intensified biological pump in the first decade, acting through regionally coherent physical (MLD-mediated mixing) and ecological (size-structure shift) pathways.

The coupling between NCP and phytoplankton biomass during the second decade aligns with expectations, as SIC remained relatively stable, reflecting post-melt stage conditions. During this stage, NCP was mainly driven by nutrient availability and primary production, consistent with mechanisms observed in the open ocean. Variability in NPP during this decade was largely influenced by increased phytoplankton biomass, resulting from the intensified inflow of nutrient-rich waters from the subarctic ocean—a process commonly referred to as “Atlantification” or “Pacification”. This process strengthened the coupling between NCP and Chla (Fig. 3e and Supplementary Fig. 13c). The lack of a significant correlation between the e-ratio and Chla is not unexpected, because the simultaneous enhancement of NCP and NPP, due to increased phytoplankton biomass, suppressed variability in their ratio (i.e., the e-ratio), thereby weakening the correlation.

We further applied a structural equation model (Piecewise SEM) to quantify the causal pathways linking sea ice decline, mixed-layer processes, and biological production (Fig. 4). Results showed that SIC exerted the strongest influence on both NCP and the e-ratio (Fig. 3d and Supplementary Fig. 13a, b), confirming sea ice decline as the dominant driver of biological pump trends during the rapid melt period of 2003–2012. During this period, SIC decline affected NCP and the e-ratio through both direct and indirect pathways. These effects were mediated by biogeochemical processes, including MLD deepening and nutrient entrainment (path coefficient: β = 0.818, P < 0.001, Fig. 4a). In contrast, during 2013–2022, the model revealed a driver chain dominated by biological adaptation, in which sea ice-driven processes weakened, nutrient supply and biomass formed a positive feedback loop, and biomass (Chla) directly promoted NCP (path coefficient: β = 0.715, P < 0.05, Fig. 4b). These results imply that the primary drivers of the biological pump shifted over the past two decades, transitioning from sea ice dynamics to phytoplankton biomass as the dominant factor.

To further elucidate changes in NCP from satellite observations, we extended our analysis over the 251-year period using an ESM. Simulations were performed for a historical period (1850–2014, BHIST) and a high-emission future scenario (2015–2100, BSSP585) using the most recent release of the Community Earth System Model version 2 (CESM2) (see Methods). Although CESM2 has been used previously to simulate Arctic sea ice, snow, clouds and precipitation, and sea surface pCO2, here we extend its application to NCP modeling. To better represent plankton diversity and size structure in NCP simulations, we incorporated the Size-based Plankton Ecological TRAits (SPECTRA) model into the Marine Biogeochemistry Library (MARBL) framework (Supplementary Table 4).

Model results corroborate satellite observations, showing NCP regime shifts closely tied to SIC variability (Fig. 5a). In the early industrial period (1850s–1960s), NCP gradually increased under high and stable SIC conditions (Fig. 5b). From the 1970s to the 2020s, NCP rose sharply, coinciding with rapid sea ice decline (Fig. 5c). However, as sea ice continues to retreat, NCP is projected to decline from the 2030s, with a substantial drop between the 2080s and 2100 (Fig. 5d). Strong correlations between NCP and SIC persisted throughout the entire period. The model also showed significant long-term shifts in phytoplankton community structure (Supplementary Fig. 16a). Prior to the 1970s, phytoplankton community composition remained stable (Supplementary Fig. 16b). Between the 1970s and 2020s, large phytoplankton expanded while smaller phytoplankton declined, coinciding with higher NCP. Conversely, from the 2030s to 2100, large phytoplankton are projected to decrease, and smaller phytoplankton are expected to increase, paralleling the anticipated decline in NCP.

CESM2 simulations revealed strong spatial coherence between trends in NCP and large phytoplankton (Supplementary Fig. 17), highlighting the role of phytoplankton community structure in modulating the Arctic biological pump. For instance, regions with strong NCP increases (e.g., Barents Sea, Chukchi Sea) coincided spatially with high abundances of large phytoplankton under contemporary conditions (1981–2020). Under future projections (2021–2100), both NCP and large phytoplankton biomass are projected to decrease (Supplementary Fig. 17c, f), suggesting a potential “miniaturization” of the phytoplankton community under sustained low-SIC conditions. This shift toward smaller cells favors remineralization in the euphotic zone, reducing export efficiency and decoupling biomass from carbon sequestration. Together, satellite observations and model simulations indicate that the biological pump in the AO is highly responsive to sea ice change. Notably, the model projects another regime shift starting in the 2030s, when SIC drops below 40%. In some regions, satellite-derived NCP and the e-ratio data suggest that this shift may already be underway (Supplementary Fig. 10).

Methods
Data sources and processing
Satellite observations, field measurements, and numerical models were used to construct a spatially and temporally continuous NCP dataset, following the procedure outlined in Supplementary Fig. 1. Underway sea surface pCO2 observations in the AO were compiled from a variety of international sources, including the Surface Ocean CO2 Atlas (SOCAT v.2022; http://www.socat.info), the Global Ocean Data Analysis Project (GLODAPv2.2022; https://glodap.info/), the Japan Agency for Marine-Earth Science and Technology (JAMSTEC; http://www.godac.jamstec.go.jp/darwin/e), the Chinese National Arctic and Antarctic Data Center (Chinese NAADC; https://datacenter.chinare.org.cn/) and the United States Geological Survey database (USGS; https://pubs.er.usgs.gov). Only summer (June–September) pCO2 measurements from 2003 to 2022 were used to calibrate and validate the retrieval model for estimating sea surface pCO2 (Supplementary Fig. 2). Ultimately, the dataset includes over 456,000 sea surface pCO2 measurements collected during 193 cruises, accompanied by corresponding SST and sea surface salinity (SSS) records (Supplementary Fig. 2b). Note that the SOCAT v.2022 and GLODAPv2.2022 datasets provide the fugacity of CO2 (fCO2) rather than pCO2. To ensure comparability, data were standardized to pCO2. Values originally reported as fCO2 were converted to pCO2 based on in situ SST measurements using Eq. (1):

pCO2 = fCO2 × 1.00436 × 10−4.669 × 10−5 × SST (1)

The conversion from fCO2 to pCO2 introduces a negligible error, as the observed difference is smaller than the ±2 μatm measurement precision. Each sampling point was assigned to one of eleven subregions and four functional topographies for regional analysis (Supplementary Dataset 1). Subregional boundary maps and geolocation data (Fig. 1a and Supplementary Dataset 1) are provided to ensure reproducibility of regional assignments. In this study, MODIS/Aqua summer climatological SST data were used to divide the Nordic Sea into the Greenland Sea and the Norwegian Sea. For the sake of brevity, this will not be reiterated below.

MODIS/Aqua L3 SST data (R2019.0), including daily, 8-day, and monthly composites at a 4-km resolution, were acquired from the NASA Ocean Color website (https://oceancolor.gsfc.nasa.gov/l3/). Satellite-observed SST was compared with concurrent underway SST observations within the same spatial grid. In the AO, the root mean square error (RMSE) between satellite-observed and in situ SST measurements was 1.08 °C, with a coefficient of determination (R2) of 0.87 (see Supplementary Fig. 3a, b). A linear correction was applied to the satellite-observed SST data (Supplementary Fig. 3a) to reduce bias in subsequent analyses.

We obtained MODIS/Aqua L2 remote sensing reflectance (Rrs) products (R2022.0) from the NASA Ocean Color website (https://oceancolor.gsfc.nasa.gov/) to develop a localized Chla retrieval algorithm, calibrated using the in situ Chla dataset provided by Lewis and Arrigo. The locally calibrated Chla algorithm yielded an RMSE of 0.21 mg m–3, a MAPE of 39.6%, a Slope of 0.72, and an R2 value of 0.65 (see Supplementary Fig. 4). This algorithm was then applied to daily MODIS/Aqua imagery to estimate daily Chla concentrations across the AO. Pixels flagged for potential errors—such as atmospheric correction failure (bit 0), land (bit 1), high or saturated radiance (bit 2), stray light contamination (bit 8), and cloud/ice presence (bit 9)—were masked in daily imagery to maintain data quality. Daily data were then aggregated to generate 8-day and monthly composites.

Daily and monthly data of SSS and MLD, provided at a spatial resolution of 0.083° × 0.083° (~8 km), were compiled from the Copernicus Marine Environmental Monitoring Service (CMEMS) global ocean reanalysis product (https://data.marine.copernicus.eu/product/GLOBAL_MULTIYEAR_PHY_001_030/). The reanalysis dataset provides multi-decadal, three-dimensional fields of ocean physical properties, generated by integrating satellite and in situ observations into a numerical ocean modeling framework, thus ensuring consistent and homogeneous representations of the ocean state over time. Prior to incorporating the SSS product into the retrieval model, reanalysis SSS values were compared with field measurements, demonstrating strong consistency. Although overestimation was noted at low SSS levels (<25 psu) (Supplementary Fig. 3c), the correlation (R2 = 0.79, P < 0.001) between measured and reanalysis SSS was significant, with most residuals below 2 psu and an RMSE of 1.91 psu (see Supplementary Fig. 3d). As waters with SSS below 25 psu account for a minor proportion of the AO, a simple regression equation was applied to regions with higher salinity (>25 psu) for correction to minimize error propagation. Near-weekly xCO2 data at a 0.05 sine-of-latitude resolution were sourced from the National Oceanic and Atmospheric Administration (NOAA) Earth System Research Laboratory (ESRL) Greenhouse Gas Marine Boundary Layer (MBL) Reference (http://www.esrl.noaa.gov/gmd/ccgg/mbl/).

Model construction of pCO2 and assessment (Step 1)
Monthly SST, xCO2, Chla, SSS, and MLD were utilized to reconstruct long-term sea surface pCO2 in the AO. SST was used to capture the thermodynamic effects, xCO2 was employed to reflect the forcing of sea surface pCO2 by rising atmospheric CO2 concentrations, Chla served as a proxy for the biological effects on sea surface pCO2, and SSS and MLD were used to characterize the physical mechanisms of horizontal and vertical seawater mixing, respectively. All environmental variables were resampled to a 4-km × 4-km spatial resolution to ensure consistency.

The matchup process between field-measured variables and satellite observation-based products was conducted using consistent spatial and temporal criteria. Specifically, 8-day composites of MODIS/Aqua Chla and SST were applied for matchup with field measurements, thereby reducing cloud contamination and increasing the number of valid pairings. For each in situ sea surface pCO2 measurement, valid satellite data within a 3 × 3-pixel window were extracted and averaged to reduce sensor-related and algorithm-derived noise. Matchups were retained only when at least five valid pixels were present within the window, and the variance among these pixels did not exceed 10%. To reduce the influence of skewness arising from their multi-order variability, Chla and MLD were log-transformed (base 10) prior to analysis. We compiled 90,623 in situ–satellite matchup samples of sea surface pCO2 and associated environmental variables, which were utilized to train and validate the retrieval scheme.

Using these data, an XGBoost machine learning model reconstructed monthly pCO2 distributions across the AO for 2003–2022. During model calibration, 60% of the 90,623 matchup pairs (57,092 matchup pairs) were randomly selected to construct the training datasets for the XGBoost model, 30% of the original dataset (24,468 matchup pairs) was designated as the corresponding test dataset, and the remaining 10% (9063 matchup pairs) were reserved for independent validation of the model’s performance. The final XGBoost algorithm consisted of 1820 decision trees with a maximum depth of nine levels. Each tree was trained using 90% of the input features and 80% of the training samples, selected through random subsampling. A learning rate of 0.2 was applied to each sequential tree, and L1 regularization was used to mitigate overfitting. Multiple metrics were used to comprehensively assess the XGBoost model, including RMSE, R², mean bias (MB), mean ratio (MR), mean absolute percentage error (MAPE), and mean relative error (MRE).

Overall, the XGBoost model training results yielded an average R2 of 0.98 and an RMSE of 5 μatm compared to the observations. The MB was lower than –0.02 μatm, the MR was near 1, the MAPE was less than 1.0%, and the MRE was below –0.04%. The statistical metrics indicate that the XGBoost regression model reliably reproduced the variability of sea surface pCO2 across the AO. Its predictive skill was further assessed with an independent validation dataset (Supplementary Table 1). For the two cruises conducted in the Barents and Chukchi Seas, the XGBoost model successfully reproduced the variation of pCO2 on the shelves (Supplementary Fig. 5), with R2 of 0.95 and 0.91, RMSE of 13.87 μatm and 16.53 μatm, MB of 0.37 μatm and 1.76 μatm, MR of 1.00 and 1.01, MAPE of 2.53% and 3.03%, and MRE of 0.16% and 0.62%, respectively. For the dataset comprising 25 independent cruises, the XGBoost model achieved R2 > 0.95, RMSE < 9 μatm, an MB of –0.03 μatm, an MR near 1, MAPE less than 1.50%, and an MRE of –0.02% (Supplementary Table 1). These results confirm that the satellite-based XGBoost model reliably reconstructs sea surface pCO2 variability.

Calculation air-sea CO2 flux (Step 2)
Using satellite-derived pCO2, 10 m wind speed (U10), xCO2, sea level pressure (Psl), ice, SSS, and SST, we estimated long-term air-sea CO2 fluxes in the AO from 2003 to 2022. The net air-sea CO2 flux (FCO2) was calculated as follows:

FCO2 = KCO2 × Ks × (1 − ice) × ΔpCO2 (2)

KCO2 = 0.251 × U10^2 × (Sc / 660)^−0.5 (3)

ΔpCO2 = pCO2sea − pCO2air (4)

pCO2air = xCO2 × (Psl − Pw) (5)

where KCO2 (m d−1) in Eq. (2) and (3) represents the gas transfer velocity of CO2, estimated from the 10 m wind speed (U10, m s−1) and the Schmidt number (Sc) following Wanninkhof et al. Monthly mean U10 values were taken from the NCEP-DOE Reanalysis 2 dataset (https://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis2.html). The Schmidt number (Sc) was calculated from SST (°C) using the coefficients reported by Wanninkhof et al. Ks in Eq. (2) represents the CO2 solubility coefficient, which depends on SST and SSS. The term ice denotes the fraction of the air-water interface covered by sea ice (%, between 0 and 1), extracted from the NOAA/NSIDC Climate Data Record of Passive Microwave Sea Ice Concentration product (G02202). ΔpCO2 is the difference in partial pressure between oceanic and atmospheric CO2 (equation (4)), where pCO2sea and pCO2air represent sea surface and atmospheric pCO2, respectively. In Eq. (5), Psl and Pw are sea level pressure and water vapor pressure, respectively. Monthly mean Psl values were extracted from a reanalysis dataset (ECMWF-ERA5, https://cds.climate.copernicus.eu/) over the study period, and Pw was calculated from Psl and SST. The pCO2air values were computed from monthly averages of the xCO2 (data adopted from the above) and atmospheric pressure following Dickson et al. Similarly, all parameters, including pCO2, U10, and sea ice, were mapped onto a 4 km × 4 km grid and averaged to monthly means to ensure consistency with the gridded sea surface pCO2 as described above.

Calculation of net community production (Step 3-6)
Satellite observations have been widely used to derive NCP via various empirical algorithms, enabling further investigations of biological pump variability in both regional and global oceans. However, the inherent limitations of these empirical approaches constrain their application in the rapidly changing AO, as they either assume a constant e-ratio or rely on statistical regressions between NCP and NPP. To overcome these limitations, we developed a semi-analytical retrieval algorithm to characterize the biological pump in the changing AO and applied it to generate a satellite-based summer NCP dataset spanning nearly two decades (Supplementary Fig. 1).

NCP was estimated based on observed temporal variations of nutrients or DIC in the upper mixed layer. In this study, we employed DIC drawdown in the upper mixed layer to estimate NCP. Monthly DIC values were calculated from satellite-derived pCO2 (pCO2sat), total alkalinity (TAsat) and other parameters (Psl, SSTsat, SSSsat) using the CO2SYS program version 3.2 (Step 4). The pCO2sat values were retrieved from satellite observations utilizing the XGBoost model described above, while SSTsat represents the MODIS/Aqua L3 SST product. The SSSsat, provided by the CMEMS, is the reanalysis product simulated from satellite observations and field measurements. The TAsat in surface seawater was estimated from SSSsat using an empirical statistical relationship calibrated with discrete field measurements of TA and salinity from the upper mixed layer (depth < 30 m), obtained from the GLODAPv2.2022 dataset (Step 3, Supplementary Fig. 1c). Only samples collected in the AO during summer (June–September) from 2003 to 2022 were utilized to establish an empirical algorithm between SSS and TA. This linear relationship (Eq. (6), R2 = 0.91, N = 1664, P < 0.001) was employed to derive the monthly TAsat dataset for the AO over the long-term series from the SSS product.

TA = (44.27 ± 0.34) × SSS + (767.39 ± 10.17) (6)

Although the temporal scales of the two datasets (monthly DICsat dataset vs daily DICinsitu dataset) are not aligned, the DICsat generated by the CO2SYS program effectively captured the variation of field-measured DICinsitu, exhibiting a robust linear relationship (Slope = 0.81, R2 = 0.52, P < 0.001, Supplementary Fig. 6a). Additionally, a reference dataset (DICcal) was produced using field-measured pCO2, SSS, SST, Psl, and TA processed in the CO2SYS program to validate the satellite-based DICsat. The comparison yielded strong agreement, with a Slope of 0.88 and an R2 of 0.78 (Supplementary Fig. 6b). These results confirm that the satellite observation-based DICsat successfully captured the DIC variations and could be reliably used to estimate NCP in the AO.

A one-dimensional biogeochemical mass balance model was applied to the mixed layer to quantify the DIC budget and disentangle the contributions of physical mixing, air-sea CO2 fluxes, and biological activities (e.g., photosynthesis and respiration) (Step 5). Therefore, based on the conservation of mass principle, the temporal change in DIC mass balance (ΔDIC) can be expressed as Eq. (7):

ΔDIC = ΔDICa-s + ΔDICbio + ΔDICmix (7)

ΔDIC = DICt+1 − DICt (8)

where ΔDIC represents the total change in DIC inventory within the surface mixed layer, calculated as DIC concentration at time step t+1 (one month in this study) minus that at time step t (Eq. (8)). ΔDICa-s denotes the contribution from air-sea CO2 flux, ΔDICbio represents biologically driven DIC changes, and ΔDICmix accounts for physical mixing processes, including horizontal advection, vertical entrainment/detrainment, and turbulent diffusion.

To address horizontal advection, we applied a non-zero end-member salinity normalization approach (Eq. (9)), where salinity was used as a conservative tracer to correct the dilution effect of sea ice meltwater and river runoff on DIC:

sDIC = DIC − (DICs=0 × S / S0) + DICs=0 (9)

where the end member salinity, DIC, and TA values for meltwater/river runoff are 0, 400 μmol kg−1, and 104 μmol kg−1. Using the previously described monthly SSSsat product, we calculated the grid-cell mean salinity over 2003–2022 as a reference value (S0), which helps to smooth short-term advective variability.

In the AO, vertical turbulent exchange across the base of the mixed layer is typically weak during summer, driven by strong stratification. Reported vertical diffusion coefficients (Kz) range from 10−8 to 10−6 m2 s−1, exerting a much smaller influence on DIC freshwater inputs (i.e., sea ice meltwater and river runoff). Consequently, their contribution to DIC variability is negligible. Regarding the entrainment/detrainment process, its accurate quantification relies heavily on spatiotemporally continuous observations of DIC below the mixed layer (DICbelow). However, due to limitations in the penetration depth and inversion accuracy of remote sensing techniques, the vertical fine structure of deep DIC—such as the dynamic variations in DICbelow with depth and time—is challenging to precisely retrieve via remote sensing. To assess the potential impact of the entrainment/detrainment process, we used in situ DIC profile datasets (e.g., GLODAPv2.2022) to calculate the vertical DIC gradient (dDIC/dz) and derive DICbelow (i.e., DIC below the MLD). Results demonstrate that the entrainment/detrainment process contributes approximately 8.57% to the sDIC estimate. Additionally, MLD remains relatively stable during summer in the AO, such that vertical entrainment events are episodic. Therefore, the impact of vertical mixing on DIC was relatively minor. The temporal changes in sDIC (ΔsDIC) primarily reflect biological activity and air-sea CO2 flux (Eq. (10)):

ΔsDIC = ΔDICa-s + ΔDICbio (10)

Surface seawater density, ρ (kg m−3) was derived from SST and SSS data using the equation of state and used to convert DIC from units of μmol kg−1 to mmol m−3. Changes in DIC attributable to air-sea CO2 exchange were calculated according to Eq. (11):

ΔDICa-s = FCO2 × Δt / MLD (11)

where FCO2 represents the net air-sea CO2 flux (mmol C m−2 d−1), Δt is a timestep (one month in this study), and MLD data were obtained from CMEMS. Biologically driven DIC change (ΔDICbio) was calculated by correcting the air-sea CO2 flux effect (ΔDICa-s) based on conservative mixing (Eq. (12)). Finally, NCP (mmol C m−2 d−1, i.e., ΔDICbio) was estimated as Eq. (13):

ΔDICbio = ΔsDIC − FCO2 × Δt / MLD (12)

NCP = −ΔDICbio × MLD / Δt (13)

Using DICsat to compute ΔDIC in Eq. (7), NCPsat can subsequently be estimated from satellite observation-based products through Eq. (8–13). The reliability of satellite-derived NCP estimates (NCPsat) was assessed by comparison with ship-based NCP derived from O2/Ar measurements collected during five cruises across the western AO, central and eastern North American Arctic, and Subarctic (Supplementary Fig. 7). In high-latitude regions, ocean color data acquisition is challenged by sea ice, persistent cloud cover, and high solar zenith angles, resulting in substantial data loss after quality control. To ensure consistency, we evaluated the NCPsat and NCPinsitu measurements in space and time. Supplementary Fig. 7 shows that NCPsat and NCPinsitu exhibit strong consistency in overall spatial patterns and temporal variations. Statistical comparison across NCPsat and NCPinsitu indicates that the mean deviation of NCP ranges from −7.83 to 2.29 mmol C m−2 d−1, which falls within the expected accuracy range for NCP estimates. Although this agreement holds for most observations, several instantaneous, localized productivity peaks—especially in coastal waters—captured by shipboard measurements are underestimated in the monthly averaged NCPsat data. This underestimation mainly results from mismatches in spatial and temporal resolution. Satellite-derived values represent horizontal averages over 4 km × 4 km pixels and monthly timescales, while shipboard measurements are instantaneous, small-scale “snapshots”. Residual analysis (right panel of Supplementary Fig. 7) indicates that large residuals, though occasionally present, are rare and contribute negligibly to the total variance. The dominant pattern shows a tight clustering of residuals around zero, indicating strong overall consistency between the datasets. In summary, despite limitations in resolving short-lived, localized productivity events, NCPsat reliably captures the dominant spatiotemporal variability of NCP in the AO and is suitable for assessing large-scale, long-term trends.

Calculation of net primary productivity
We calculated NPP using five distinct models to evaluate biological pump efficiency and its ranges (Supplementary Fig. 8). The Vertically Generalized Production Model (VGPM), which is based on chlorophyll, available light, and photosynthetic efficiency, links chlorophyll to temperature-dependent growth rates via a polynomial function of SST. Although it lacks explicit spectral, temporal, or vertical resolution, VGPM remains widely applied, with studies highlighting the critical need for accurate estimation of photosynthetic efficiency. Accordingly, Morel et al. modified the photosynthetic efficiency in VGPM to an exponential function of SST based on the Eppley curve, yielding a reasonable estimate of the photosynthetic efficiency (e.g., Eppley-VGPM). Pabi et al. developed a regionally calibrated NPP algorithm, the Arctic Ocean Productivity Model (AOPM), which integrates satellite remote sensing data of Chla, SST, PAR, and sea ice cover to account for the changing AO. In contrast, the Carbon-based Productivity Model (CbPM) uses satellite-derived particulate organic carbon concentration to serve as the primary metric instead of chlorophyll. The CAFE model explicitly considers the absorption by Colored, Dissolved, and Detrital material (CDM) to reduce overestimations of NPP in coastal waters. Given the reliance on the variable Chla in multiple models (VGPM, Eppley-VGPM, and AOPM), we developed a locally optimized satellite-based Chla remote sensing retrieval algorithm (Supplementary Fig. 4) tailored for the AO to improve the accuracy of NPP estimation.

To ensure robust trend analysis, we examined long-term temporal trends across all five NPP datasets (Supplementary Fig. 11a). VGPM, Eppley-VGPM, and AOPM showed broadly similar trajectories over the two decades. In contrast, CbPM continued to rise significantly after 2012, while the CAFE trend weakened, reflecting inherent algorithmic sensitivities to environmental drivers or input data sources. Importantly, spatial distribution patterns were highly consistent across all products (Supplementary Fig. 8). The e-ratio (NCP/NPP) remained remarkably stable both spatially and temporally (Supplementary Figs. 8 and 11a), demonstrating robustness against NPP model choice and supporting its use in assessing export efficiency. To further minimize inter-model discrepancies in trend analysis, we used the mean across the five NPP products (Supplementary Fig. 9). We focused on NPP and NCP rates (mmol C m−2 d−1) rather than total integrated carbon (Tg C) to disentangle ecological dynamics from the confounding influence of expanding ice-free areas. Previous studies have indicated that total Arctic NPP (Tg C) increases primarily reflect ice-free area expansion rather than changes in NPP rate. Our multi-product trend analysis confirms the absence of significant trends in NPP rate during the first decade (Supplementary Fig. 11a), consistent with stable Chla biomass proxies.

Uncertainty analysis
Uncertainties in pCO2 and FCO2 products were quantified following the methodological approach of Landschützer et al., as implemented in later studies. For each grid cell, the total uncertainty in estimated pCO2 was composed of three sources: the direct pCO2 measurement uncertainty (Uobs) from the databases, gridding uncertainty (Ugrid), and mapping uncertainty (Umap). In accordance with Sharp et al., Uobs was assigned based on SOCAT fCO2 accuracy flags (A-D), with flags “A” and “B” corresponding to ±2 µatm, and flags “C” and “D” reflecting an accuracy of 5 µatm. Only data from the SOCAT database corresponding to Flags “A” and “B” were employed for constructing the retrieval model. The uncertainties of pCO2 measurements from other databases remain unknown, and their sample sizes (2.52%) are limited. To ensure a conservative approach, Uobs was set at a maximum of 5 μatm. The uncertainty derived from comparing pCO2 products generated by the XGBoost model with observations from the pCO2 databases represents Umap, with a predictive error of 8.63 μatm in the AO. Ugrid, representing the uncertainty introduced when gridding pCO2 observations from the databases onto 4-km monthly mesh maps, was calculated at 1.82 µatm, based on the estimates from the method of Roobaert et al. for the AO. Assuming independence among these sources, the overall uncertainty of the gridded pCO2 product (UpCO2) was determined using Monte Carlo error propagation (Supplementary Table 2).

In estimating FCO2, uncertainties arise from four independent components in equations (2) to (5) used to calculate FCO2: (1) the XGBoost-derived sea surface pCO2 field constraining the air–sea gradient, (2) the gas exchange transfer velocity, (3) the solubility coefficient of carbon dioxide, and (4) the choice of the sea-ice product. Assuming the independence of these four sources, the overall FCO2 uncertainty (θFCO2) was derived through root-sum-of-squares propagation (Eq. (14)), following the strategy utilized by Landschützer et al.:

θFCO2 = 2 * sqrt(θΔpCO2^2 + θkCO2^2 + θks^2 + θice^2) (14)

the uncertainty of the transfer velocity, θkCO2, primarily associated with the wind product, was estimated to be approximately 20%. The contribution of the gas solubility coefficient term (θks) to our uncertainty assessment is minimal (0.2%), as suggested by Weiss et al. Assessments by Peng et al. indicate that the uncertainty associated with sea ice products is roughly 5%. To ensure a conservative estimate, we doubled this value and adopted 10% as the sea ice uncertainty. The NOAA MBL reference product provides an estimate of xCO2 uncertainty, derived by averaging uncertainties across each sine latitude step in the AO. As reported by Lan et al., the average uncertainty for xCO2 is 0.32 ppm. Additional uncertainties for other variables are detailed in Supplementary Table 2. Subsequently, we employed Eqs. (4) and (5) to calculate the uncertainties for pCO2air and ΔpCO2. Because NCP uncertainty is derived through the propagation of errors in DIC, we first evaluated the sources of DIC uncertainty. Given the uncertainties of relevant parameters (i.e., salinity-based TA, pCO2; Supplementary Table 2), the propagated uncertainties of DIC were estimated using an add-on for the CO2SYS program in MATLAB 2021. The propagated uncertainty incorporates combined standard uncertainties from input variables as well as the equilibrium constants. TA uncertainty was assessed prior to the CO2SYS calculations and subsequent error propagation. The uncertainty of DIC was estimated at ~37.52 μmol kg−1, based on uncertainties in the TA calculation (41.44 μmol kg−1; see Eq. (6)), pCO2 (0.31 μatm), and outputs from CO2SYS. We also split the AO into four subregions (Inflow shelf, Interior shelf, Outflow shelf, and Basin), and developed empirical algorithms for each region using salinity as a proxy to estimate the surface TA concentrations, yielding an average uncertainty of 38.38 μmol kg−1, with a total uncertainty in the calculated DIC of 34.95 μmol kg−1. This result is slightly lower than the DIC uncertainty calculated from the salinity-TA relationship (Eq. (6)).

We quantified NCP uncertainties (equations (7) to (13)) by propagating parameter error through Monte Carlo simulations. For each error term, 100,000 simulations were run using Gaussian distributions, and the 95% confidence interval was used to represent uncertainty bounds for each data point in the time series. Combined uncertainty from individual sources of error as well as NCP using the Monte Carlo simulation approach is shown in Supplementary Table 2. Satellite and reanalysis data contributed the most to the total uncertainty. The gas transfer coefficient uncertainty θkCO2 (20%), represents the largest contributor to FCO2 uncertainty, which propagates further into NCP calculation. The resulting uncertainty for simulated NCP based on the Monte Carlo approach is 46.85%. All reported uncertainties apply to the pan-Arctic scale and may be larger at the subregional scale. The uncertainty associated with in situ NCP measurements was negligible (0.2%) relative to other error terms, and therefore excluded from the final uncertainty analysis. Although the uncertainty in NCP approaches 50%, these estimates still discern the direction of NCP measurements observed during most cruises. Moreover, this uncertainty does not significantly affect the spatial distribution or long-term trend of NCP, but only its magnitude.

NCP simulation based on Community Earth System Model version 2 (CESM2)
We conducted the historical scenario experiment (1850−2014) and the high-emissions future scenario (SSP5-8.5) experiment (2015−2100) using CESM2 configured with MARBL-SPECTRA to assess longer-term changes in NCP and its relationship with SIC across the pan-Arctic region. CESM2 integrates multiple components, including the CAM6 atmosphere, CLM5 land, MOSART river runoff, WW3 wave, CISM2 land-ice, CICE5 sea-ice, POP2 ocean, and MARBL ocean biogeochemistry models. The Marine Biogeochemistry Library (MARBL) is the core oceanic biogeochemical component of the CESM2. Its modular design allows for flexible configuration of plankton functional types. In this study, we implemented the Size-based Plankton Ecological TRAits (SPECTRA) community model within MARBL. This configuration, referred to as MARBL-SPECTRA, represents an advancement over the default MARBL setup, including nine phytoplankton and six zooplankton functional types with equivalent spherical diameter (ESD) ranging from 0.47–300 μm and 0.02–20 mm, respectively. Compared to the default configuration (three types of phytoplankton and one type of zooplankton), MARBL-SPECTRA significantly expands the diversity of functional types and size classes of plankton. This enables a more complex representation of the marine food web and improves our understanding of key processes such as plankton community dynamics, which were previously difficult to resolve.

MARBL-SPECTRA supports multiple ecosystem configurations of varying complexity, each incorporating different quantities and properties of functional groups, thereby demonstrating the flexibility of MARBL’s ecosystem configuration. We first used spectra-config, available at https://github.com/marbl-ecosys/spectra-config, to generate configuration files and initial conditions for MARBL-SPECTRA. MARBL-SPECTRA was then integrated into CESM2, available at https://github.com/ESCOMP/CESM, to simulate long-term biogeochemical variables from 1850 to 2100. The experiments were conducted with BHIST and BSSP585 configurations at f19_g16 resolution, which uses a finite-volume 1.9°×2.5° mesh for the atmospheric and land components and a gx1v6 polar-shifted grid (~1°) for the oceanic and sea-ice modules. All remaining CESM2 parameters were kept at their default settings. From the CESM2 outputs, the rate of change in DIC due to biological activity was used to assess NCP variability in the AO. Sea ice fraction was also extracted to examine the variation of SIC. After modifying the phytoplankton community module in the CESM2 model, we extracted biomass data for nine phytoplankton species (Supplementary Table 4), grouped them into small and large categories based on ESD, and calculated their relative proportions in the AO. The small group (ESD < 20 μm) included the two smallest mixed phytoplankton (mp1, mp2), picoplankton (pp), and diazotroph (diaz), whereas the large group (ESD ≥ 20 μm) comprised three diatom size classes (diat1, diat2, diat3) and the two largest mixed phytoplankton types (mp3, mp4). The extracted data, including NCP, sea ice fraction, and phytoplankton community structure, were subsequently used to analyse the relationships among these variables in the AO.

Statistical analysis
We applied Sen’s slope (derived from the Mann-Kendall test) of the long-term summer means of NCP, NPP, and the e-ratio to estimate their rates of change, with the corresponding P values used to assess the significance of these trends. Kendall’s tau-b correlation coefficient was used to assess correlations between time and NCP, NPP, e-ratio, SIC, and OWD, with trends deemed significant at P < 0.05. All Mann-Kendall trend analyses were conducted in MATLAB 2021. Additionally, Pearson’s correlation analyses (Pearson’s r) were performed to evaluate relationships between summer mean NCP and e-ratio and variables including SIC, OWD, Chla, and MLD. Statistical significance of these correlation coefficients was determined using P values (e.g., P < 0.05). Additionally, to investigate the causal relationship between key drivers and two metrics of the Arctic biological carbon pump—NCP and carbon export efficiency (e-ratio), we categorized potential drivers into physical variables (SIC, Wind, MLD, OWD, PAR, and SST) and biological variables (Nutrient, Chla, and NPP). All variables were derived from the dataset described in the Methods section. Nutrient data were obtained through modal decomposition of nitrate, phosphate, and silicate concentrations, with the original data sourced from the CMEMS product. Informed by prior knowledge, we employed the “Piecewise SEM” R package to construct a piecewise structural equation model, quantifying the pathways through which these drivers affect NCP and the e-ratio. The final model was selected based on the goodness-of-fit statistic (Fisher’s C) and the corresponding P value.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Using satellite-based datasets of NCP and a numerical model, this study offers a more systematic understanding of the Arctic biological pump, particularly in terms of its spatial variability and multi-decadal changes—surpassing previous insights that were based on limited in situ measurements. The authors report that the biological pump intensity, measured as the e-ratio, exhibits distinct trends over the past two decades. Specifically, they observed a rapid increase during the first decade (2003–2012), followed by no significant trend in the second decade (2013–2022), independent of NPP products. They attribute this shift to a change in primary drivers from sea ice loss and its cascading physical effects in the first decade to biological factors of phytoplankton structure in the second. The numerical model provides consistent evidence supporting this transition and helps elucidate the underlying mechanisms. The authors clearly describe their methodology and conduct careful data processing, including uncertainty analysis and data correction. Overall, this is a well-written manuscript that provides a novel perspective on the responses of biological pump to the rapidly changing Arctic Ocean.\n\nMy main concern is the insufficient comprehensive on explaining the primary drivers behind the decadal trends in the e-ratio. First, the reported shift in the e-ratio requires rigorous statistical analysis to assess its significance. Do similar shifts occur in the key driving variables, and do their temporal changes align with the e-ratio trend? Currently, the authors rely primarily on time series correlations to propose the underlying mechanisms, which may not be sufficient to establish causality. While the spatial variations in the e-ratio are well described, there is a lack of corresponding analysis on the spatial variability of MLD and phytoplankton community structure. I recommend including spatial correlation analyses between the e-ratio and MLD and phytoplankton structure to enhance the robustness and credibility of the proposed mechanisms.\n\nAnother concern lies in the presented NPP trends, which appear inconsistent with previous studies (e.g., ref 43 in this manuscript) that reported a significant increase in NPP associated with sea ice loss, followed by a shift to phytoplankton biomass-driven increases in recent decade. These earlier findings indicate a similar temporal transition in NPP and its dominant drivers, coinciding with the period identified in this study. How do the authors explain this discrepancy? I note that the authors use an average of five distinct NPP products but do not present a clear increasing trend in NPP during the first decade. Is this trend consistently absent across all products, or do some exhibit significant increases while others do not? Clarifying this would help reconcile the differences with previous literature and strengthen the overall interpretation.\n\nOther comments:\n1.Line 62: How is the role of “vertically migrating zooplankton” reflected in this study? It would be helpful to clarify whether and how their contribution to the biological pump is explicitly or implicitly represented in the analysis.\n2.Lines 201–203: Is the deepening of the mixed layer more pronounced in regions where the e-ratio shows a faster increase? A spatial comparison or statistical relationship would help clarify this link.\n3.Lines 207–209: Can a change in MLD of just ~2 meters realistically lead to such a large change in NCP? Please elaborate on the underlying mechanism driving the correlation between MLD and NCP.\n4.Lines 217 and 250: It would be beneficial to provide an explanation of how larger phytoplankton contribute to increased NCP. A literature-based discussion would enhance the interpretation.\n5.I note a near-perfect correlation between NCP and SIC based on the CESM2 model in Fig.4b-d. What is the governing equation for DIC in CESM2? I wonder whether the parameterization of biological processes may be overly sensitive to SIC changes. Additionally, what is the modeled relationship between NPP and SIC?\n6.Equation (7): In considering changes in the DIC inventory, how is horizontal or vertical advection accounted for over the monthly time scale? Please clarify whether the advection effect has been ignored or approximated.\n7.Figure 7: What explains the missing portions of the satellite-based NCP data? High values of underway O2/Ar-based NCP appear underrepresented. A monthly mean-state versus snapshot comparison may not be statistically appropriate. A scatter plot of all underway measurements versus satellite estimates could be overly noisy—would it be more informative to use daily satellite-derived NCP for the comparison?\n8.Figure 9: It is recommended to include the long-term trends of all five NPP products in this figure to allow readers to assess consistency and variability across different datasets.", "reviewer_2": "Review of Sea Ice Loss Leads to Regime Shifts in the Arctic Biological Pump’ by Wu, Hu et al.\n\nThis manuscript uses a mix of observations and model output to quantify trends in Arctic Ocean productivity and carbon export this century (for obs) and from 1850 to 2100 (for models). In general, I think the manuscript uses an appropriate but technically demanding method to quantify variability in biological productivity and carbon export over the last ~2 decades. The authors convincingly link the variability to changes in sea ice (however see my detailed comments for questions about the relative role of mixing, light and sea ice).\n\nI recommend publication after revision.\n\nThe introduction, lines 36-78, sets up the problem well. Sea ice has declined but there is disagreement over the impact on the biological carbon pump, because inferences are drawn from sparse observations. Satellite-based or gridded products can help resolve this. The authors then briefly describe their method, which is explained in fine detail in the methods section, before moving on to results.\n\nResults (bio pump strength and efficiency): The main point of this manuscript is the trends in NCP and export and their relationship to changes in sea ice. I found the description of regional variability (lines 102-130) less interesting. If the authors want to free up some space to discuss trends more, then these 2 paragraphs could mostly be deleted.\n\nResults (changes in the biological pump trends): First, I think ‘changes’ and ‘trends’ in this sub-heading are redundant, just use one or the other, not both. Importantly, the units of the trends are wrong in the manuscript text and Figure 1 (panels a, b, c). NCP and NPP have units of mmolC m^-2 year^-1, so the trends should have units of mmolC m^-2 year^-2. The units for the e ratio trend are correct in the ms text and Figure 2 panel d (year^-1), but missing in panels e and f. As above, I found the description of the regional variability not very interesting.\n\nResults (Biological pump shifting linked to sea ice decline): I’d suggest replacing ‘shifting’ with ‘trends’ in this sub-heading. This section is a good description of trends and their potential drivers. It’s also a part of the ms where regional variability in sea ice is described (eg lines 175-181). So, contrary to what I said above, it seems regional variability in the biological pump is relevant, but the two are not linked. If the authors want to discuss regional variability in NCP/e-ratio and sea-ice, they should connect them better.\n\nLinks between SIC, mixing and export, lines 198-202: Here, the authors describe changes in SIC, MLD and export. They describe a transition from a melt-water dominated regime to a mixed-layer dominated regime. These terms don’t have much meaning for me. I view them as related: melt water influences stratification and therefore MLD. Maybe what the authors mean is that mixing from storms starts to become more important than stratification driven by meltwater. Also please add some explanation about how a ‘mixed layer dominated regime’ (use a more informative term perhaps) enhances export. Same at line 208 where it is stated that reduced upper ocean stability enhances export. Just add a few words explaining how.\n\nLines 210-212: It might be just a grammar issue, but it appears that the authors are linking deeper light penetration with enhanced vertical mixing of nutrients. They are different processes, so please clarify.\n\nLine 256: Delete ‘dynamically’. It doesn’t really mean anything in this context.\n\nLine 261: I would re-title this section ‘conclusions’ rather than ‘discussion’.\n\nFigures: In general the figures are well done and informative, but I have a few suggestions.\n\nFigure 1: In the caption text, sometimes the letter for the figure panel comes before the text describing it, sometimes after. Please be consistent (this comment applies to other figures too). On line 295, I’m not sure what ‘back’ is there for.\n\nFigure 2: The units for the trends are wrong, see comments above.\n\nFigure 3: Panels d-g are described as showing ‘correlations’. It’s true that the figures contain text that quantifies the correlations, but I don’t agree that the plots themselves are showing correlations. To do that it would be better to plot the properties against each other, like a scatter plot with NPP on the y axis and SIC on the x axis, not a time series of each variable. Figure 4 more correctly describes the plots as time series, not correlations.\n\nMETHODS\nThe methods are important and well done, but I will comment on one detail. In general, despite its complexity, I like the approach of trying to use all available data to properly model the carbon system.\n\nLine 525 mentions delta-DIC_mix but I feel like the contribution to DIC from vertical mixing is not well explained. Similarly, around line 535, the authors explain how they account for the contribution to DIC from horizontal mixing of meltwater and river water, but vertical mixing is not mentioned. Please explain the role of vertical mixing (or maybe it has no role?).\n\n(Remarks on code availability)\nZenodo link did not work."}, "author_response": "Dear reviewers,\nThank you for offering us an opportunity to improve the quality of our submitted\nmanuscript (“Sea Ice Loss Leads to Regime Shifts in the Arctic Biological Pump”). We\nappreciated very much the reviewers’ constructive and insightful comments. In this\nrevision, we have addressed all of these comments/suggestions. We hope that the\nrevised manuscript has now addressed the reviewers’ concerns.\nWe highlighted all the revisions in red colour. Our responses to reviewers are explained\nbelow, in blue colour beneath each comment.\nResponse to Reviewer #1\nReviewer #1 (Remarks to the Author):\nUsing satellite-based datasets of NCP and a numerical model, this study offers a more\nsystematic understanding of the Arctic biological pump, particularly in terms of its\nspatial variability and multi-decadal changes—surpassing previous insights that were\nbased on limited in situ measurements. The authors report that the biological pump\nintensity, measured as the e-ratio, exhibits distinct trends over the past two decades.\nSpecifically, they observed a rapid increase during the first decade (2003–2012),\nfollowed by no significant trend in the second decade (2013–2022), independent of NPP\nproducts. They attribute this shift to a change in primary drivers from sea ice loss\nand its cascading physical effects in the first decade to biological factors of phytoplankton\nstructure in the second. The numerical model provides consistent evidence supporting\nthis transition and helps elucidate the underlying mechanisms. The authors clearly\ndescribe their methodology and conduct careful data processing, including uncertainty\nanalysis and data correction. Overall, this is a well-written manuscript that provides a\nnovel perspective on the responses of biological pump to the rapidly changing Arctic\nOcean.\nResponse: We sincerely thank the reviewer for the careful review and positive feedback\non our manuscript. All constructive comments and suggestions have been carefully\naddressed in the revised manuscript. A point-by-point response to each comment is\nprovided below.\nMy main concern is the insufficient comprehensive on explaining the primary drivers\nbehind the decadal trends in the e-ratio. First, the reported shift in the e-ratio requires\nrigorous statistical analysis to assess its significance. Do similar shifts occur in the key\ndriving variables, and do their temporal changes align with the e-ratio trend? Currently,\nthe authors rely primarily on time series correlations to propose the underlying\nmechanisms, which may not be sufficient to establish causality. While the spatial\nvariations in the e-ratio are well described, there is a lack of corresponding analysis on\nthe spatial variability of MLD and phytoplankton community structure. I recommend\nincluding spatial correlation analyses between the e-ratio and MLD and phytoplankton\nstructure to enhance the robustness and credibility of the proposed mechanisms.\nResponse: We have incorporated a spatial correlation analysis between the e-ratio and\nMLD in the revised manuscript (Figure R1). Results reveal that MLD exhibits\nsignificant positive correlations with NCP and the e-ratio across most shelf regions of\nthe Arctic Ocean, which further supports the time series analysis presented in Fig. 3g.\nFigure R1. Spatial correlation maps of MLD, NCP, and e-ratio for 2003-2012. (a) Correlation\nbetween MLD and NCP; (b) Correlation between MLD and e-ratio. Blue and red colours represent\nnegative and positive correlations, respectively, with stippling denoting statistically significant\nrelationships (P < 0.05).\nTo improve the generalization of our conclusions, we supplemented the spatial analysis\nwith outputs of the CESM2 model, which simulate changes in NCP and phytoplankton\ncommunity structure. Results show that regions with significant NCP increases (e.g.,\nthe Barents Sea and Chukchi Sea; Figure R2b) largely overlap with areas where large\nphytoplankton biomass has increased (Figure R2e). This overlap suggests that these\nregions exhibit synchronous responses to environmental disturbances during the early\nstages of climate change (initial sea ice melt phase). Under future scenarios (2021–\n2100), the upward trend in NCP weakens and the growth of large phytoplankton\nbiomass stagnates (Figures R2c,f). This pattern reflects a negative feedback on the\nbiological carbon pump, driven by a shift in ecological niches toward the dominance of\nsmall phytoplankton. Following carbon fixation, the small phytoplankton are more\nreadily decomposed within the euphotic zone. This decomposition lowers NCP and\ndecouples phytoplankton biomass from the efficiency of carbon fixation. Spatially,\ntrends in large phytoplankton biomass and NCP exhibit strong coupling in the ice-\nmargin zone and consistent gradients across distinct sea regions. Ecologically, the high\nefficiency of the diatom-dominated biological pump supports our conclusion that the\nproliferation of large phytoplankton directly drives increases in NCP. In summary, the\nexpansion of large phytoplankton in Arctic marginal seas is both an ecological response\nto climate change and a primary driver of NCP variability. This finding underscores the\npivotal role of phytoplankton community structure in regulating the efficiency of the\nArctic biological carbon pump (i.e., the e-ratio).\nFigure R2. Trends in NCP and large phytoplankton biomass from 1850 to 2100 derived from\nCESM2 output. NCP trends for (a) 1850-1980, (b) 1981-2020, and (c) projected 2021-2100. Large\nphytoplankton trends for (d) 1850-1980, (e) 1981-2020, and (f) projected 2021-2100. Colour scales\nindicate the magnitude and direction of trends, with blue denoting negative and red denoting positive\ntrends. Black dots mark statistically significance trends at the 5% level (P < 0.05).\nTo more rigorously investigate the causal relationships between environmental drivers\nand the e-ratio, we constructed a piecewise structural equation model (PiecewiseSEM;\nFig. R3) by incorporating well-established physiological and ecological pathways,\nfollowing the method described in Lefcheck et al.1.\nThe PiecewiseSEM quantifies the direct and indirect effects of key drivers (e.g., SIC,\nMLD) on NCP and the e-ratio, effectively distinguishes causal pathways from covariate\nones, and enables causal inference regarding the underlying driving mechanisms at the\nstatistical model level. As shown in Figure R3a, the model exhibits a good fit (Fisher's\nC = 66.032, P = 0.612, df = 70, AIC = 90.126) and identifies a sea ice–dominated\ndriving chain for 2003–2012 that melting of SIC prolongs open-water duration (OWD),\nintensifies wind-field disturbances, and consequently deepens the MLD. This\ndeepening of MLD, in turn, directly or indirectly enhances NCP, thereby increasing the\ne-ratio. Figure R3b presents a biological adaptation–dominated driver chain for 2013–\n2022, in which sea ice-driven processes weaken, nutrient supply and phytoplankton\nbiomass form a positive feedback loop, and phytoplankton biomass (Chla) directly\npromotes NCP. This finding supports our conclusion that the biological pump has\nshifted from a sea ice-dominated regime to a phytoplankton biomass-dominated regime.\nFigure R3. Direct and indirect effects of environmental drivers on NCP and the e-ratio, derived\nfrom a structural equation model (PiecewiseSEM). Blue boxes indicate physical factors, and\norange boxes indicate biological factors. Solid black arrows denote positive relationships, while\ndashed red arrows denote negative relationships. Numbers on the arrows show significant\nstandardized path coefficients, and asterisks indicate significance levels (*P < 0.05, **P < 0.01,\n***P < 0.001).\nAll these revisions have been incorporated into the revised manuscript.\nAnother concern lies in the presented NPP trends, which appear inconsistent with\nprevious studies (e.g., ref 43 in this manuscript) that reported a significant increase in\nNPP associated with sea ice loss, followed by a shift to phytoplankton biomass-driven\nincreases in recent decade. These earlier findings indicate a similar temporal transition\nin NPP and its dominant drivers, coinciding with the period identified in this study.\nHow do the authors explain this discrepancy? I note that the authors use an average of\nfive distinct NPP products but do not present a clear increasing trend in NPP during the\nfirst decade. Is this trend consistently absent across all products, or do some exhibit\nsignificant increases while others do not? Clarifying this would help reconcile the\ndifferences with previous literature and strengthen the overall interpretation.\nResponse: Good point. We have incorporated trends in NPP from multiple products\ninto the revised manuscript. Results show that these products exhibit consistent patterns\nover the past two decades (Supplementary Fig. 11a). The discrepancy between our\nfindings and those of previous studies stems from the use of different metrics: we\nemployed NPP rate (mmol C m−2 d−1), whereas Ardyna et al.2 used area-integrated NPP\n(Tg C). A similar trend would emerge if we applied the same area-integrated NPP metric\n(Figure R4).\nThis clarification has been incorporated into the revised manuscript (Lines 622–636).\nFigure R4. Trends of NPP in the AO. (a) NPP rate (mmol C m−2 d−1) and (b) area-integrated NPP (Tg\nC). The correlation coefficient (Kendall’s tau-b r, from Mann-Kendall test) indicate the long-term rate\nin NPP, with associated P value used to assess the significance of the trends are presented in the panels.\nDashed lines represent Sen’s slope of NPP versus time from different periods.\nOther comments:\n1.Line 62: How is the role of “vertically migrating zooplankton” reflected in this study?\nIt would be helpful to clarify whether and how their contribution to the biological pump\nis explicitly or implicitly represented in the analysis.\nResponse: We appreciate this insightful comment and concur that vertically migrating\nzooplankton play a crucial role in regulating the biological pump across global oceans.\nHowever, quantifying their contributions across the pan-Arctic over long-term time\nseries remains methodologically challenging.\nIn this study, we used NCP—defined as the difference between NPP and heterotrophic\nrespiration (HR)—as a metric to assess the strength of the biological pump. This metric\nencompasses all processes influencing carbon retention and transport, including those\nmediated by vertically migrating zooplankton.\nFor instance, carbon transport by vertically migrating zooplankton—via feeding on\nsurface phytoplankton and subsequent excretion or mortality at depth—is indirectly\nreflected in NCP, as it modifies the community-level balance between production (i.e.,\nNPP) and respiration (i.e., HR). Therefore, although the contribution of vertically\nmigrating zooplankton is not explicitly quantified, it is implicitly embedded in the\necosystem carbon balance represented by NCP. Notably, while NCP integrates multiple\nprocesses, it does not allow for the isolation of specific contributions from individual\nmechanisms—e.g., vertical migration of zooplankton or vertical mixing of dissolved\norganic carbon (DOC). Disentangling the individual contributions of vertically\nmigrating zooplankton and other processes (e.g., vertical mixing of DOC) is essential\nfor advancing our understanding of the regulatory mechanisms governing the biological\npump—and represents a key objective of our future research.\nThis clarification has been incorporated into the revised manuscript (Lines 71–75;\nLines 276–278).\n2.Lines 201–203: Is the deepening of the mixed layer more pronounced in regions\nwhere the e-ratio shows a faster increase? A spatial comparison or statistical\nrelationship would help clarify this link.\nResponse: Good point. Following your suggestion, we have added the spatial\ncomparison and the statistical relationship between MLD and the e-ratio in the revised\nmanuscript (see Figure R1).\nThe revised paragraph in the Results section (“Physical–biological coupling driving\nshifts in the biological pump”) now reads as follows:\n“Significant correlations between SIC or OWD and MLD (Fig. 3f; Supplementary Figs.\n12e and 13d), along with a positive correlation between NCP or e-ratio and MLD\nduring the active melting period of 2003–2012 (Fig. 3g; Supplementary Figs. 13e and\n14), support the hypothesis that reduced upper ocean stability, reflected by MLD\ndeepening under declining SIC, was the key physical driver of enhanced carbon export\n(higher NCP and higher e-ratio).”\n3.Lines 207–209: Can a change in MLD of just ~2 meters realistically lead to such a\nlarge change in NCP? Please elaborate on the underlying mechanism driving the\ncorrelation between MLD and NCP.\nResponse: We appreciate this constructive suggestion, which has assisted us in\nclarifying the potential impacts of MLD deepening on NCP and the e-ratio. A mixed\nlayer deepening of ~2 m represents a relative increase of approximately 10–40% or\nmore in the Arctic Ocean. This is because MLD is typically shallow in summer, ranging\nfrom 5 to 30 m3-5.\nAs discussed in the manuscript, such MLD deepening can substantially affect NCP via\nmultiple mechanisms. As shown in Supplementary Fig. 13, surface waters in the Arctic\nare typically nutrient-depleted, as denoted by the light-colored regions. Despite\nabundant solar radiation in summer, phytoplankton growth is strongly limited by\nnutrient availability. Environmental perturbations, such as increased wind stress\nassociated with sea ice retreat, can deepen MLD and entrain nutrient-rich subsurface\nwaters (darker red regions), thereby enhancing nutrient supply and stimulating NCP.\nUnder nutrient-limited conditions, phytoplankton communities exhibit non-linear\nresponses to nutrient enrichment, and even modest nutrient inputs can trigger\nsubstantial increases in productivity. Tremblay et al.6 reported that such nutrient\nenrichment can double NCP, while Mills et al.7 demonstrated that nitrate (NO −)\n3\nadditions increase Chla concentrations by 1.4- to 6.8-fold relative to controls.\nCorresponding increases in POC further confirm that NCP is nutrient-limited and\nhighly responsive to changes in nutrient availability.\nMoreover, phytoplankton photophysiology is highly sensitive to variations in light\navailability. Although sea ice retreat increases solar radiation by reducing surface\nalbedo, it can expose surface phytoplankton to excessive light—inducing\nphotoinhibition and reducing photosynthetic efficiency8. MLD deepening redistributes\nphytoplankton to deeper, lower-light layers—alleviating light stress, improving\nphotosynthetic efficiency, and ultimately enhancing NCP. Notably, Kellogg et al.9.\nfound that nutrient availability influences phytoplankton cell size, with higher nutrient\nlevels favoring the proliferation of larger phytoplankton taxa. Nutrient input associated\nwith MLD deepening specifically favors the growth of large, fast-growing\nphytoplankton (e.g., diatoms)—taxa that contribute disproportionately to NCP. As the\nphytoplankton community shifts toward larger taxa (Supplementary Fig. 14), NCP\nincreases substantially.\nTo further investigate the linkage between MLD and biological pump efficiency, we\nhave added a spatial correlation analysis between MLD and e-ratio (Figure R1) and\nexamined their causal relationship using a structural equation model (Figure R3). These\nadditions provide a mechanistic explanation for the observed MLD-NCP relationship.\nAdditional discussion of this linkage has been incorporated into the revised manuscript\n(Lines 183–204).\n4.Lines 217 and 250: It would be beneficial to provide an explanation of how larger\nphytoplankton contribute to increased NCP. A literature-based discussion would\nenhance the interpretation.\nResponse: Great suggestion. We have added further discussion to enhance the\ninterpretation in Lines 195–201 of the revised manuscript:\n“Larger phytoplankton enhance carbon export because their faster sinking accelerates\nPOC transfer to depth10. Rapid sinking shortens particle residence time in the euphotic\nzone and limits remineralization by heterotrophs11. As a result, a larger fraction of\nprimary production is exported as sinking particles12,13 (i.e., the e-ratio; Fig. 2d),\ncontributing to higher NCP. The link between phytoplankton size structure and carbon\nexport explains the rise in NCP despite the absence of a significant upward trend in\nNPP during this period2,14 (Supplementary Fig. 9d).”\n5.I note a near-perfect correlation between NCP and SIC based on the CESM2 model\nin Fig.4b-d. What is the governing equation for DIC in CESM2? I wonder whether the\nparameterization of biological processes may be overly sensitive to SIC changes.\nAdditionally, what is the modeled relationship between NPP and SIC?\nResponse: In CESM2, the evolution of DIC is governed by the prognostic tracer\nequation implemented in the MARBL15 (Marine Biogeochemistry Library;\nhttps://marbl.readthedocs.io/en/latest/dev-guide/introduction.html) framework:\n∂χ\n+ ∇∙(𝑢𝜒)−∇∙(𝐾 ∙∇𝜒) = Β (𝑥)\n∂t 𝜒\nHere, 𝜒 denotes DIC, the left-hand side of the equation accounts for temporal changes,\nadvection, and diffusion (calculated by the ocean general circulation model, OGCM),\nwhile the right-hand side Β𝜒 (𝑥) represents the net sources and sinks of DIC\n(calculated by MARBL). For DIC, Β𝜒 (𝑥) includes biological processes, air-sea CO\n2\nexchange, and chemical buffering. MARBL calculates these terms based on model state\nvariables (e.g., phytoplankton biomass, nutrient concentrations, and temperature).\nAlthough SIC is not explicitly included in the equation, it indirectly influences the\nbiological processes in Β𝜒 (𝑥) by modulating key environmental factors. For\nexample, a decrease in SIC increases light availability under ice-covered regions,\nthereby enhancing phytoplankton photosynthesis and biological DIC uptake; in parallel,\nsea ice melting alters the MLD, affecting the vertical transport of nutrients (e.g., nitrate\nand silicate) and consequently modulating phytoplankton growth and the net DIC\nbudget. Thus, the strong correlation between NCP and SIC observed in this study\nprimarily reflects the indirect effect of SIC on the biological budget of DIC via\nenvironmental regulation. Meanwhile, NCP serves as an indicator of the rate at which\nDIC is consumed by biological activity in the upper ocean.\nTo address the concern regarding whether biological processes are overly sensitive to\nSIC variability, we compared two configurations of the MARBL module:\n(1) Official MARBL configuration: includes three phytoplankton functional groups—\ndiatoms, diazotrophs, and “small” (pico/nano) phytoplankton—with simplified\nparameterization to simulate core ecological functions;\n(2) Modified MARBL (MARBL-SPECTRA): includes nine phytoplankton functional\ngroups, expanding particle size spectra and nutrient acquisition strategies, and refining\nphysiological parameters to better capture ecological diversity.\nFigure R5. Time series of summer NCP and SIC simulated by the MARBL module based on the\nofficial CESM2.\nDespite fundamental differences in community complexity and physiological\nparameterization, both configurations simulate a strong response of NCP to changes in\nSIC (see Figure R5). This finding suggests that the strong NCP-SIC correlation is not\nan artifact of a specific biological parameterization or overly sensitive biological\nparameters, but rather reflects a robust feature of the model’s physical response to sea\nice variability.\nModeled NPP also exhibits a significant and persistent correlation with SIC, with a\ntrend that closely mirrors the NCP-SIC relationship. However, the NPP-SIC correlation\nweakens after 2030 (though it remains statistically significant), whereas the NCP-SIC\ncorrelation remains consistently strong. This divergence highlights CESM2’s capacity\nto capture the complexity and nonlinearity of sea ice-ecosystem interactions, while also\nsupporting the existence of distinct biogeochemical feedback pathways across different\nindicators. While NPP represents the “production side”, NCP integrates the full\n‘production-respiration-export’ chain, thereby enhancing the mechanistic interpretation\nand credibility of the model outputs. The consistent correlations between NCP-SIC and\nNPP-SIC suggest that the model captures fundamental biogeochemical responses to sea\nice variability, rather than artifacts of parameterization schemes.\nFigure R6. Time series of summer NPP and SIC simulated by CESM2. (a) Temporal variability of\nNPP and SIC from 1850 to 2100. Trends in NPP and SIC during the (b) historical (1850–1970),\n(c) modern (1970–2030), and (d) future periods (2030–2100), respectively. In (a), blue, yellow, and red\nshading indicate distinct trend regimes of NCP and SIC. Thin and bold lines denote annual values and\nfive-year running means of NPP and SIC, respectively. The upper-right corner of each panel shows long-\nterm changes in NPP and SIC based on Sen’s slope and the Mann-Kendall significance test (P value).\nPearson’s correlation coefficient (r) and corresponding P value for the NPP-SIC relationship are shown\nin the lower-left corner of the panel.\n6.Equation (7): In considering changes in the DIC inventory, how is horizontal or\nvertical advection accounted for over the monthly time scale? Please clarify whether\nthe advection effect has been ignored or approximated.\nResponse: Good point. In our model, both horizontal and vertical advection are\nexplicitly represented within the physical mixing term (ΔDICₘᵢₓ) in equation (7).\nSpecifically, we apply a non-zero endmember salinity normalization method16-18\n(equation (9)) to correct for DIC dilution or enrichment caused by physical mixing.\nThis approach effectively isolates errors induced by physical mixing—including\nadvection—by separating them from the influences of biological activity and air–sea\nCO exchange. Moreover, at a monthly temporal resolution (Δt = 1 month), high-\n2\nfrequency advective variability (e.g., mesoscale eddies) is smoothed out, whereas large-\nscale advection is partially captured through the climatological reference salinity field\n(S ). The strong agreement between satellite-derived NCP estimates and in-situ O /Ar-\n0 2\nbased measurements (Supplementary Fig. 7) further validates the effectiveness of our\nsalinity normalization method (equation (9)) in minimizing physical mixing-related\nuncertainties, thereby reinforcing the robustness of our monthly NCP estimates in the\nAO.\nThis clarification has been incorporated into the revised manuscript (Lines 536–546).\n7.Figure 7: What explains the missing portions of the satellite-based NCP data? High\nvalues of underway O2/Ar-based NCP appear underrepresented. A monthly mean-state\nversus snapshot comparison may not be statistically appropriate. A scatter plot of all\nunderway measurements versus satellite estimates could be overly noisy—would it be\nmore informative to use daily satellite-derived NCP for the comparison?\nResponse: We acknowledge the existence of data gaps in the satellite-derived NCP\ndataset. These gaps are primarily attributed to extensive sea ice coverage and frequent\ncloudiness in high-latitude regions19, which impede the retrieval of reliable ocean color\nproducts, such as Chla. Additionally, large solar zenith angles in high-latitude regions\nfurther compromise the reliability of satellite observations.\nTo ensure data quality, we applied strict preprocessing quality control measures (see\nMethods), excluding pixels affected by clouds, sea ice, aerosols, or large solar zenith\nangles. While these procedures improve the reliability of the retained data, they\ninevitably reduce data availability.\nSatellite-based observations also underestimate the high NCP values detected by ship-\nbased O /Ar measurements. This discrepancy is primarily driven by differences in the\n2\nspatial and temporal resolutions of the two methods. Ship-based observations provide\nhigh-resolution, localized “snapshots” along cruise transects, capturing transient, small-\nscale productivity hotspots—such as algal blooms and upwelling zones—with high\nsensitivity. In contrast, satellite-based NCP products rely on coarser spatial grids (4 km\nin this study) and represent area-averaged values, which inherently smooth extreme\nproductivity peaks.\nFurthermore, satellite inversion algorithms—particularly those based on Chla—can\nexhibit non-linear saturation under high-biomass conditions, which reduces their\nsensitivity to extreme NCP values. Additionally, shipboard measurements capture near-\ninstantaneous variability (minutes to hours), whereas satellite products, averaged over\nlonger timescales (monthly composites), obscure brief productivity spikes driven by\nrapid bloom development or mixing events. Consequently, even with spatial alignment,\ntemporal averaging in satellite products constrains the detection of peak NCP intensities\nobserved in-situ.\nWe concur that using daily satellite-derived NCP for the comparison would be more\ninformative. However, persistent cloud cover and sea ice in the Arctic render valid daily\nsatellite observations extremely scarce (Figure R6). This scarcity results in few matched\ndata points within the intended comparison windows, thereby reducing both the\nfrequency of matching and statistical robustness. Therefore, to ensure sufficient sample\nsize for meaningful validation, we used monthly-averaged satellite NCP products,\nwhich provide broader spatial and temporal coverage for this study.\nFigure R7. Comparison between satellite-derived daily NCP and underway O /Ar-based NCP\n2\nmeasurements from multiple cruises.\nThe primary objective of our satellite NCP algorithm is to leverage the extensive spatial\ncoverage and long-term observational continuity of remote sensing to characterize\nspatial patterns and seasonal-to-interannual variability in NCP. The underestimation of\nhigh ship-based NCP values may introduce regional uncertainties. However, this\nlimitation is likely to have minimal impact on our long-term, pan-Arctic scale analysis.\nFurthermore, a comprehensive uncertainty analysis confirms that satellite-based\nobservations capture key spatial and seasonal patterns in NCP across the Arctic Ocean\n(Figure R8; Figure S7 in the revised manuscript).\nThis additional analysis is now detailed in our manuscript (Lines 580–599).\nFigure R8. Comparison of satellite-derived NCP with underway O2/Ar-based NCP measurements\nfrom multiple cruises. Sample distribution collected by the RV Xuelong during two Chinese National\nArctic Research Expedition (CHINARE) cruises conducted from (a) July 24 to September 4, 2016, and\n(c) July 29 to September 8, 2018, covering the western AO between 65°N–85°N and 137°W–180°W.\nO2/Ar-based NCP measurements collected from the CCGS Amundsen during ArcticNet expeditions in\n(e) 2015 (10 July–20 August), (g) 2018 (24 July–8 September), and (i) 2019 (5 July–15 August), covering\nthe central and eastern North American Arctic and Subarctic. Monthly satellite-derived NCP estimates\nwere paired with contemporaneous underway O2/Ar-based NCP measurements. Spatial distributions of\nO2/Ar-derived NCP along cruise tracks are shown in the subpanels. Values for the root mean square error\n(RMSE), mean absolute error (MAE), and mean bias (MB) are reported. The histogram on the right of\neach panel displays the residual distribution between satellite-derived and underway measurements\nacross the cruises. Gray curves indicate the fitted probability density function.\n8.Figure 9: It is recommended to include the long-term trends of all five NPP products\nin this figure to allow readers to assess consistency and variability across different\ndatasets.\nResponse: We thank the reviewer for the insightful comments and constructive\nsuggestions, which have substantially enhanced the quality of our manuscript. The\nlong-term temporal trends of all five NPP products have been added to the revised\nSupplementary Fig. 11a.\nWe have incorporated additional details in the revised manuscript (Lines 622–636).\nReviewer #2 (Remarks to the Author):\nReview of Sea Ice Loss Leads to Regime Shifts in the Arctic Biological Pump’ by Wu,\nHu et al.\nThis manuscript uses a mix of observations and model output to quantify trends in\nArctic Ocean productivity and carbon export this century (for obs) and from 1850 to\n2100 (for models). In general, I think the manuscript uses an appropriate but technically\ndemanding method to quantify variability in biological productivity and carbon export\nover the last ~2 decades. The authors convincingly link the variability to changes in sea\nice (however see my detailed comments for questions about the relative role of mixing,\nlight and sea ice).\nI recommend publication after revision.\nResponse: We thank the reviewer for the supportive comments on our manuscript. All\nconstructive suggestions have been carefully considered and addressed in the revised\nmanuscript. Please find our point-by-point responses to your comments below.\nThe introduction, lines 36-78, sets up the problem well. Sea ice has declined but there\nis disagreement over the impact on the biological carbon pump, because inferences are\ndrawn from sparse observations. Satellite-based or gridded products can help resolve\nthis. The authors then briefly describe their method, which is explained in fine detail in\nthe methods section, before moving on to results.\nResponse: Thanks for your positive comments.\nResults (bio pump strength and efficiency): The main point of this manuscript is the\ntrends in NCP and export and their relationship to changes in sea ice. I found the\ndescription of regional variability (lines 102-130) less interesting. If the authors want\nto free up some space to discuss trends more, then these 2 paragraphs could mostly be\ndeleted.\nResponse: Good suggestion. In the revised manuscript, we have condensed and merged\nthe two paragraphs into a single paragraph (Lines 106–122) to succinctly describe the\nregional variability in NCP and the e-ratio.\nResults (changes in the biological pump trends): First, I think ‘changes’ and ‘trends’ in\nthis sub-heading are redundant, just use one or the other, not both. Importantly, the units\nof the trends are wrong in the manuscript text and Figure 2 (panels a, b, c). NCP and\nNPP have units of mmolC m^-2 year^-1, so the trends should have units of mmolC m^-\n2 year^-2. The units for the e ratio trend are correct in the ms text and Figure 2 panel d\n(year^-1), but missing in panels e and f. As above, I found the description of the regional\nvariability not very interesting.\nResponse: We appreciate these constructive points. The subtitle has been revised to\n“Trends in the Biological Pump” in the revised manuscript. Units in the figures have\nbeen updated, and the description of regional variability has been condensed (Lines\n125–142).\nResults (Biological pump shifting linked to sea ice decline): I’d suggest replacing\n‘shifting’ with ‘trends’ in this sub-heading. This section is a good description of trends\nand their potential drivers. It’s also a part of the ms where regional variability in sea ice\nis described (eg lines 175-181). So, contrary to what I said above, it seems regional\nvariability in the biological pump is relevant, but the two are not linked. If the authors\nwant to discuss regional variability in NCP/e-ratio and sea-ice, they should connect\nthem better.\nResponse: We appreciate your valuable suggestions. The subtitle has been revised to\n“Biological Pump Trends Linked to Sea Ice Decline” in the updated manuscript.\nAdditionally, we have added several sentences to better link regional variability in sea\nice concentration (SIC) and the biological pump in this paragraph (Lines 147–156),\nas follows:\n“Between 2003 and 2022, Sea Ice Concentration (SIC) declined significantly while\nopen water duration (OWD) increased (Fig. 3a and Supplementary Fig. 12a),\ncoinciding with NCP changes. From 2003 to 2012, SIC declined rapidly, and NCP\nincreased significantly, particularly in interior and outflow shelves such as the Kara,\nLaptev, and northern Barents Seas, as well as the Beaufort Sea (Fig. 3b and\nSupplementary Fig. 12b). In contrast, from 2013 to 2022, the rate of SIC loss slowed,\nleaving concentrations relatively stable. In this period, SIC decreased slightly—but not\nsignificantly—in the eastern interior shelves (East Siberian, Laptev, and Kara Seas),\nwhile a slight increase was observed in the western shelves (Beaufort Sea) (Fig. 3c and\nSupplementary Fig. 12c). Correspondingly, the biological pump remained stable (Fig.\n3d).”\nLinks between SIC, mixing and export, lines 198-202: Here, the authors describe\nchanges in SIC, MLD and export. They describe a transition from a melt-water\ndominated regime to a mixed-layer dominated regime. These terms don’t have much\nmeaning for me. I view them as related: melt water influences stratification and\ntherefore MLD. Maybe what the authors mean is that mixing from storms starts to\nbecome more important than stratification driven by meltwater. Also please add some\nexplanation about how a ‘mixed layer dominated regime’ (use a more informative term\nperhaps) enhances export. Same at line 208 where it is stated that reduced upper ocean\nstability enhances export. Just add a few words explaining how.\nResponse: Great point. We have revised the term “melt-water dominated regime” to\n“meltwater-stratified regime”, following Vappen et al.20. This term describes conditions\nwith high sea-ice concentration, where low-salinity meltwater accumulates at the sea\nsurface during the melting season, generating strong salinity-driven stratification.\nThese conditions result in a shallow mixed layer (i.e., reduced MLD), limited vertical\nexchange in the upper ocean, and a weakened biological carbon pump. The term\n“mixed-layer dominated regime” has also been revised to “mixed-layer-homogenized\nregime”, which refers to periods of low SIC characterized by the expansion of open\nwater and the intensification of wind-driven turbulence. These conditions disrupt\nmeltwater-induced stratification, leading to a deepened mixed layer (increased MLD)\nand substantially enhanced vertical exchange in the upper ocean.\nThese revisions to terminology and their associated definitions have been clarified in\nthe revised manuscript (Lines 174–177).\nLines 210-212: It might be just a grammar issue, but it appears that the authors are\nlinking deeper light penetration with enhanced vertical mixing of nutrients. They are\ndifferent processes, so please clarify.\nResponse: Yes, this misunderstanding was caused by a grammatical error. We intended\nto state that reduced SIC significantly increases sunlight penetration into the water\ncolumn and promotes nutrient supply from deeper layers through vertical mixing and\nstorm-induced upwelling on shallow shelves. This has been clarified and revised in the\nmanuscript (Lines 187–195). The updated text now reads as follows:\n“For example, MLD deepening entrains nutrient-rich bottom waters, alleviating\nnutrient limitations in strongly stratified surface layers (see Supplementary Fig. 15).\nNutrient-deficient phytoplankton exhibit disproportionate productivity gains from small\nnutrient inputs6,7. Reduced SIC also markedly increases light penetration into the water\ncolumn21. Deeper mixing redistributes phytoplankton out of the surface light-inhibition\nzone and optimizes photosynthetic efficiency under enhanced light conditions.\nEnhanced nutrient supply9, combined with optimized photosynthesis, supports higher\nprimary production in marginal shelves2,14,22 and drives a shift in phytoplankton\ncommunities toward large-celled taxa such as diatoms23-26.”\nLine 256: Delete ‘dynamically’. It doesn’t really mean anything in this context.\nResponse: We have deleted ‘dynamically’ as suggested.\nLine 261: I would re-title this section ‘conclusions’ rather than ‘discussion’.\nResponse: We agree and have updated.\nFigures: In general the figures are well done and informative, but I have a few\nsuggestions.\nResponse: Thank you for the positive comment.\nFigure 1: In the caption text, sometimes the letter for the figure panel comes before the\ntext describing it, sometimes after. Please be consistent (this comment applies to other\nfigures too). On line 295, I’m not sure what ‘back’ is there for.\nResponse: All figure captions have been revised. The word “back” was a typo error and\nshould have referred to panel “a”. This error has been corrected in the revised\nmanuscript. We appreciate your careful review.\nFigure 2: The units for the trends are wrong, see comments above.\nResponse: The units were updated to mol C m−2 year−2, as reflected in the revised\nmanuscript (Line 311).\nFigure 3: Panels d-g are described as showing ‘correlations’. It’s true that the figures\ncontain text that quantifies the correlations, but I don’t agree that the plots themselves\nare showing correlations. To do that it would be better to plot the properties against\neach other, like a scatter plot with NPP on the y axis and SIC on the x axis, not a time\nseries of each variable. Figure 4 more correctly describes the plots as time series, not\ncorrelations.\nResponse: We agree that it is incorrect to describe panels d–g in Figure 3 as\n“correlations”. The primary purpose of this figure is to illustrate the long-term temporal\nvariability of different environmental variables in the Arctic Ocean and their\ncoincidence with NCP variations. The figure caption has been revised to clarify this.\nHowever, the correlation coefficients are retained to indicate correlations between NCP\nand the environmental variables. Additionally, we analyzed these relationships using\nscatter plots, as suggested, in the revised manuscript (Figure R9; Supplementary Fig.\n13; Line 1324).\nFigure R9. Relationships between NCP, the e-ratio and environmental factors over two periods\n(2003–2012 and 2013–2022). (a) SIC vs NCP, (b) SIC vs e-ratio, (c) Chla vs NCP, (d) SIC vs MLD, and\n(e) MLD vs NCP. Points are colored to indicate year (2003–2022). Solid lines show linear regressions\nfor each period, with correlation coefficients (Pearson’s r) and P values reported for each relationship.\nMETHODS\nThe methods are important and well done, but I will comment on one detail. In general,\ndespite its complexity, I like the approach of trying to use all available data to properly\nmodel the carbon system.\nResponse: We sincerely appreciate your positive assessment of our methodological\napproach.\nLine 525 mentions delta-DIC_mix but I feel like the contribution to DIC from vertical\nmixing is not well explained. Similarly, around line 535, the authors explain how they\naccount for the contribution to DIC from horizontal mixing of meltwater and river water,\nbut vertical mixing is not mentioned. Please explain the role of vertical mixing (or\nmaybe it has no role?).\nResponse: Good suggestion. Vertical mixing could affect the DIC inventory through\nvertical turbulent diffusion27 and entrainment/detrainment. In the AO, vertical turbulent\nexchange across the base of the mixed layer is typically very weak during summer due\nto strong stratification. Reported vertical diffusion coefficients28 (𝐾 ) range from 10−8\n𝑍\nto 10−6 m2 s−1, with their effect on DIC being several orders of magnitude smaller than\nthat of horizontal freshwater inputs (i.e., sea ice meltwater and river runoff).\nConsequently, their contribution to DIC variability is negligible.\nRegarding the entrainment/detrainment process, its accurate calculation is highly\ndependent on spatiotemporally continuous data of DIC below the mixed layer\n(DIC )29. However, due to limitations in the penetration depth and inversion\nbelow\naccuracy of remote sensing techniques, the vertical fine structure of deep DIC—such\nas the dynamic variations in DIC below the MLD and time—is challenging to precisely\nretrieve via remote sensing. To assess the potential impact of the\nentrainment/detrainment process, we used in-situ DIC profile datasets (e.g., GLODAP)\nto calculate the vertical DIC gradient (dDIC/dz) and derive DIC (i.e., DIC below\nbelow\nthe MLD)29. Results demonstrate that the entrainment/detrainment process contributes\napproximately 8.57% to the sDIC estimate. Additionally, MLD remains relatively\nstable during summer, such that vertical entrainment events are episodic. Therefore, the\nimpact from vertical mixing on DIC was relatively minor.\nThis has been clarified in the revised manuscript (Lines 536–563).\nReviewer #2 (Remarks on code availability):\nZenodo link did not work.\nResponse: The relevant scripts and results from this study have been uploaded to our\nZenodo repository: https://doi.org/10.5281/zenodo.14020827.\nWe sincerely thank the reviewer for the careful review and supportive comments. These\ninsightful suggestions have been invaluable in enhancing the quality of our manuscript.\nReference\n1. Lefcheck, J. S. & Freckleton, R. piecewiseSEM: Piecewise structural equation modelling inr for\necology, evolution, and systematics. Methods Ecol. Evol. 7, 573-579 (2015).\n2. Ardyna, M. & Arrigo, K. R. Phytoplankton dynamics in a changing Arctic Ocean. Nat. Clim.\nChang. 10, 892-903 (2020).\n3. Brenner, S. et al. Wind-Driven Motions of the Ocean Surface Mixed Layer in the Western Arctic.\nJ. Phys. Oceanogr. 53, 1787-1804 (2023).\n4. Hordoir, R. et al. Changes in Arctic Stratification and Mixed Layer Depth Cycle: A Modeling\nAnalysis. J. Geophys. Res.-Oceans. 127, e2021JC017270 (2022).\n5. Allende, S., Fichefet, T., Goosse, H. & Treguier, A. M. On the ability of OMIP models to\nsimulate the ocean mixed layer depth and its seasonal cycle in the Arctic Ocean. Ocean Model.\n184, 102226 (2023).\n6. Tremblay, J. É. et al. Vertical stability and the annual dynamics of nutrients and chlorophyll\nfluorescence in the coastal, southeast Beaufort Sea. J. Geophys. Res.-Oceans. 113, C07S90\n(2008).\n7. Mills, M. M. et al. Nitrogen Limitation of the Summer Phytoplankton and Heterotrophic\nProkaryote Communities in the Chukchi Sea. Front. Mar. Sci. 5, 362 (2018).\n8. Lewis, K. M. et al. Photoacclimation of Arctic Ocean phytoplankton to shifting light and\nnutrient limitation. Limnol. Oceanogr. 64, 284-301 (2018).\n9. Kellogg, D. R. & Levin, P. A. Nutrient availability as an arbiter of cell size. Trends Cell Biol.\n32, 908-919 (2022).\n10. Miettinen, T. P. et al. Cell size, density, and nutrient dependency of unicellular algal\ngravitational sinking velocities. Sci Adv 10, eadn8356 (2024).\n11. Guidi, L. et al. A new look at ocean carbon remineralization for estimating deepwater\nsequestration. Glob. Biogeochem. Cycle 29, 1044-1059 (2015).\n12. Dunne, J. P., Armstrong, R. A., Gnanadesikan, A. & Sarmiento, J. L. Empirical and mechanistic\nmodels for the particle export ratio. Glob. Biogeochem. Cycle 19, GB4026 (2005).\n13. Dunne, J. P., Sarmiento, J. L. & Gnanadesikan, A. A synthesis of global particle export from the\nsurface ocean and cycling through the ocean interior and on the seafloor. Glob. Biogeochem.\nCycle 21, GB4006 (2007).\n14. Lewis, K. M., van Dijken, G. L. & Arrigo, K. R. Changes in phytoplankton concentration now\ndrive increased Arctic Ocean primary production. Science 369, 198-202 (2020).\n15. Long, M. C. et al. Simulations with the Marine Biogeochemistry Library (MARBL). J. Adv.\nModel. Earth Syst. 13, e2021MS002647 (2021).\n16. Ouyang, Z. et al. Seasonal water mass evolution and non-redfield dynamics enhance CO uptake\n2\nin the Chukchi Sea. J. Geophys. Res.-Oceans. 127, e2021JC018326 (2022).\n17. Ouyang, Z. et al. Sea-ice loss amplifies summertime decadal CO increase in the western Arctic\n18. Joshi, A. P. & Warrior, H. V. Comprehending the role of different mechanisms and drivers\naffecting the sea-surface pCO and the air-sea CO fluxes in the Bay of Bengal: A modeling\n2 2\nstudy. Mar. Chem. 243, 104120 (2022).\n19. Cesana, G. V. et al. The correlation between Arctic sea ice, cloud phase and radiation using A-\nTrain satellites. Atmos Chem Phys. 24, 7899-7909 (2024).\n20. von Appen, W.-J. et al. Sea-ice derived meltwater stratification slows the biological carbon\npump: Results from continuous observations. Nat. Commun. 12, 7309 (2021).\n21. Sumata, H. et al. Regime shift in Arctic Ocean sea ice thickness. Nature 615, 443-449 (2023).\n22. Hill, V., Light, B., Steele, M. & Sybrandy, A. L. Contrasting sea-ice algae blooms in a changing\nArctic documented by autonomous drifting buoys. J. Geophys. Res.-Oceans. 127,\ne2021JC017848 (2022).\n23. Blais, M. et al. Contrasting interannual changes in phytoplankton productivity and community\nstructure in the coastal Canadian Arctic Ocean. Limnol. Oceanogr. 62, 2480-2497 (2017).\n24. Toullec, J. et al. Processes controlling aggregate formation and distribution during the Arctic\nphytoplankton spring bloom in Baffin Bay. Elementa-Sci. Anthrop. 9, 1 (2021).\n25. Nadaï, G., Nöthig, E.-M., Fortier, L. & Lalande, C. Early snowmelt and sea ice breakup enhance\nalgal export in the Beaufort Sea. Prog. Oceanogr. 190, 102479 (2021).\n26. Lee, Y. et al. Influence of sea ice concentration on phytoplankton community structure in the\nChukchi and East Siberian Seas, Pacific Arctic Ocean. Deep-Sea Res. Part I-Oceanogr. Res.\nPap. 147, 54-64 (2019).\n27. Schultz, C. et al. Modeling phytoplankton blooms and inorganic carbon responses to sea-ice\nvariability in the West Antarctic Peninsula. J. Geophys. Res.-Biogeosci. 126, e2020JG006227\n(2021).\n28. Cole, S. T. & Roemer, P. A. The Transition Layer and Remnant Transition Layer of the Western\nArctic Ocean: Stratification, Vertical Diffusivity, and Pacific Summer Water Heat Fluxes. J.\nGeophys. Res.-Oceans. 129, e2023JC020059 (2024).\n29. Yang, B., Shadwick, E. H., Schultz, C. & Doney, S. C. Annual Mixed Layer Carbon Budget for\nthe West Antarctic Peninsula Continental Shelf: Insights From Year‐Round Mooring\nMeasurements. J. Geophys. Res.-Oceans. 126, e2020JC016920 (2021)."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have thoroughly addressed the concerns and points raised in my review, making thoughtful revisions and offering additional clarifications where appropriate. These changes have substantially improved the manuscript, and I am pleased to recommend it for publication.", "reviewer_2": "I have read the response to reviewers document and I am satisfied that the authors have addressed my (Reviewer 2) concerns. Thanks for such a clear explanation of changes made."}, "author_response": ""}]

Decision letter