North Atlantic subtropical mode water formation controlled by Gulf Stream fronts

ABSTRACT The North Atlantic Ocean hosts the largest volume of global subtropical mode waters (STMWs) in the world, which serve as heat, carbon and oxygen silos in the ocean interior. STMWs are formed in the Gulf Stream region where thermal fronts are pervasive and result in feedback with the atmosphere. However, their roles in STMW formation have been overlooked. Using eddy-resolving global climate simulations, we find that suppressing local frontal-scale ocean-to-atmosphere (FOA) feedback leads to STMW formation being reduced almost by half. This is because FOA feedback enlarges STMW outcropping, attributable to the mixed layer deepening associated with cumulative excessive latent heat loss due to higher wind speeds and greater air-sea humidity contrast driven by the Gulf Stream fronts. Such enhanced heat loss overshadows the stronger restratification induced by vertical eddies and turbulent heat transport, making STMW colder and heavier. With more realistic representation of FOA feedback, the eddy-present/rich coupled global climate models reproduce the observed STMWs much better than the eddy-free ones. Such improvement in STMW production cannot be achieved, even with the oceanic resolution solely refined but without coupling to the overlying atmosphere in oceanic general circulation models. Our findings highlight the need to resolve FOA feedback to ameliorate the common severe underestimation of STMW and associated heat and carbon uptakes in earth system models.

3 Atmospheric Physics (IAP) product [2]. Oceanic variables, including temperature and salinity, from all three datasets are objectively-analyzed monthly in situ observations with a horizontal resolution of 1°×1° and vertically interpolated into 10-m intervals using the Akima's scheme [3].
Considering the bias due to sparseness of in situ ocean observations before the 1970s [4], here we calculated STMW in EN4 and IAP for a climatological 30-year period from 1981 to 2010.
To evaluate the FOA feedback intensity, we used the monthly SST and air-sea net heat flux derived from the fifth-generation reanalysis product (ERA5) [5] produced by the European Centre for Medium-Range Weather Forecasts (ECMWF). It has been demonstrated that the highresolution (horizontal grid spacing of 0.25°) ECMWF ERA5 reanalysis enables the representation of oceanic frontal-scale signals along with oceanic fronts' and eddies' imprints on the overlying atmosphere [6]. The period of 1981-2010 was used for computing the climatological mean state of the FOA feedback intensity around the STMW formation region.

STMW definition
The definition of STMW used in this study aims to identify the pycnostad with a low potential vorticity (PV) constraint in contrast to the strong stratification in the pycnocline. Here, PV = − 0 , where f is the Coriolis parameter, 0 is a reference density (1024 kg m -3 ), is the potential density of seawater, and z is the vertical coordinate, positive upwards. The STMW computed from the Argo and EN4 products and CESM twin simulations is defined as a pycnostad between 26.2-26.6 kg m -3 with a PV of < 1×10 -10 m -1 s -1 and a thickness of > 100 m. For the IAP product the potential density constraint of 26.2-26.8 kg m -3 was used to encompass the low-PV water mass ( Fig. S2).
We depicted the potential density and other properties of STMW by tracking its core layer, which is defined as a local vertical minimum PV within the STMW. The core layer preserves the 4 winter mixed layer conditions and is regarded as the best record of water mass properties at the time of formation. Following ref. 1, the core layer properties are derived by volume-weighted average of the monthly values in the STMW distribution region during May-to-December when the STMW properties are not destructed apparently by outcropping.

Heat content budget in the upper ocean
The heat content budget for the upper ocean can be derived as: where h is the depth for vertical integration, T is the potential temperature, u = (uh, w) is the threedimensional velocity vector, with uh = (u, v) being its horizontal component and w being its vertical component,  is the horizontal biharmonic diffusion coefficient, 0 and are the ocean reference density and heat capacity, respectively. TD represents the heat content tendency, Qmf represents the heat transport convergence by the mean flows, Qeddyh represents the horizontal eddy heat transport convergence, Qeddyv represents the vertical eddy heat transport, Qturb represents the parameterized microscale turbulent mixing, and Qnet represents the surface net heat flux (positive into the upper ocean). The heat content budget is conducted for the upper 300 m with the corresponding variables being interpolated onto a regular 10-m vertical grid. Given the negligible effect of horizontal turbulent mixing on the upper ocean heat content [7,8], the last term is not considered in the budget analysis. The overbar denotes the monthly mean values, which are available in the model outputs. The prime denotes the deviation from the monthly mean value, 5 which can be derived by subtracting the monthly mean values from the model's diagnostic outputs.
<…> denotes a spatial average over the STMW formation region plus a time average in the winter season (ONDJFM) from year 46 to year 55.

Water mass formation framework
Walin formalism. We identified and quantified key mechanisms for the difference in STMW production with and without FOA feedback by means of water mass formation framework of Walin [9]. Before this analysis, we verified that the STMW as well as its seasonal cycle can be well defined by the core layer density (σ = 26.47 kg m -3 and σ = 26.36 kg m -3 ) along with the density windows of ∆σ = ± 0.1 kg m -3 , even without additional constraints on PV or other properties (Fig. S10). This satisfies the prerequisite for exploring mechanisms underlying STMW formation in the Walin framework, which focuses on water masses within a particular density class with no restrictions on stratification [10].
The annually integrated volume budget [11] can be written as: where V is the volume of STMW bounded by the aforementioned isopycnal surfaces σ ± ∆σ and where Qnet is the surface net heat flux as the sum of sensible heat flux, latent heat flux, and radiative heat flux, S is surface salinity, E is evaporation, P is precipitation, R is runoff, ⃗ is wind stress vector, ⃗⃗ is the unit upward vector, b is buoyancy ( = − ∆ / 0 ), is gravitational acceleration, , , are the heat capacity, thermal expansion, and haline contraction coefficient of seawater, respectively. Positive buoyancy flux implies ocean buoyancy gain (i.e., a decrease in surface density).
The water mass transformation rate TR by the air-sea fluxes is calculated by integrating B over the outcrop window bounded by isopycnal surfaces σ ± ∆σ/2. The discretized expression is given by: where ds is an area element at the sea surface. The instantaneous map of transformation can be expressed by 0 ( ∆ ) −1 ( , ∆ ) with ( , ∆ ) being the top-hat function of which is zero except within the interval of ∆σ/2 , where it has unit value [12]. Positive transformation corresponds to volume flux towards increasing density. The convergence of TR drives net 7 subduction across the surface into the ocean interior and yields the rate at which a given density class ( 1 < < 2 ) is created or destroyed via air-sea buoyancy flux (i.e., the formation rate FR): where ′ = − , ′ = − , the subscript C and F denote CTRL and FILT, respectively. The terms on the right-hand side of equation (6) represent the contributions from differences in LHF and surface outcrops between CTRL and FILT, respectively. The surface outcrops difference depends on differences in outcrop areas within the window centered at the lower and upper bounds of STMW density between CTRL and FILT.

LHF reconstruction
We reconstructed LHF from surface wind speed (U) and air-sea humidity contrast (dq) using the following bulk formulae [13]: where is the atmospheric surface density, is the latent heat of vaporization as a function of surface temperature ( ), expressed as = (2.501 − 0.00237 × ) × 10 6 , is the stability dependent turbulent exchange coefficient for latent heat, is the saturation specific humidity at , and is the atmospheric specific humidity at the lowest level. Here, the coefficients are all Based on the LHF reconstruction, the LHF difference between CTRL and FILT can be approximated as: In general, models with finer oceanic resolution are typically accompanied by higher atmospheric 9 resolution. The 15 multi-resolution simulations used in this study include 6 eddy-free, 7 eddypresent, and 2 eddy-rich simulations (Table S1). Given different density stratifications in each model, we adjusted the PV threshold and the potential density constraint to define STMWs in individual models appropriately.
The improvement of STMW simulation with increasing ocean resolution alone based on oceanic general circulation models (OGCMs) participating in Ocean Model Intercomparison Project phase 2 (OMIP-2) were used to make comparison with that in CGCMs in order to understand the role of FOA feedback in STMW production. The three OGCMs used in this study (Table S2) are forced by the Japanese 55-year atmospheric reanalysis (JRA-55) [14] in the absence of its response to oceanic feedback, covering the period from 1958 to 2018. Please refer to ref. [15] for more details of the experimental design.

FOA feedback intensity
To measure the intensity of FOA feedback in CGCMs, we regressed the wintertime (ONDJFM) monthly frontal-scale net heat flux onto the frontal-scale SST: where the regression coefficient denotes the FOA feedback intensity. The prime denotes the frontal-scale anomaly computed as the spatially high-pass filtered field achieved by removing a 5°×5° boxcar running mean. We focused on inter-model differences in the total STMW volume and its relationship with the wintertime mean FOA feedback intensity, area-averaged in the key region (29°-41°N, 75°-38°W) where mode water is produced and FOA feedback is very strong across models. We confirmed that the main results in this study are insensitive to modest changes in either the spatial extent of the region used for area averaging or the box size used for spatial filter, consistent with the fact that remains relatively stable on spatial wavelengths < 1000 km [16]. The AWI-CM-1-1 was not used in the comparison of FOA feedback intensity because of unavailability of heat flux data.

Bootstrap test
A bootstrap method [17] is used to examine whether the differences of STMW volume budget terms between CESM CTRL and FILT are statistically significant. Each term is resampled randomly to conduct 10,000 realizations of mean value for CTRL and FILT (Fig. 3a). In this resampling process, any annual mean value can be selected more than once. The standard deviation of the 10,000 inter-realizations of mean value for each simulation is then computed. If the mean value difference between CTRL and FILT is greater than the sum of those two standard deviations, then it is statistically significant above the 95% confidence level.

Taylor diagram
We used Taylor  the spatial correlation coefficient, the spatial standard deviation, and the centered pattern rootmean-square (CRMS) difference that incorporates the evaluation of similarities of spatial structure and variation between model and observations. Both the spatial standard deviation and CRMS are normalized by the observational standard deviation to facilitate direct comparison. Therefore, lying closer to the observations indicates a more realistically simulated STMW in CGCMs.

Percentage change
The percentage change in STMW thickness, volume, and formation rate when FOA feedback is suppressed are calculated as: ( − )/ × 100%. The percentage changes based on      Table S1 for further information.