The Arctic Subpolar Gyre sTate Estimate: Description and Assessment of a Data‐Constrained, Dynamically Consistent Ocean‐Sea Ice Estimate for 2002–2017

A description and assessment of the first release of the Arctic Subpolar gyre sTate Estimate (ASTE_R1), a data‐constrained ocean‐sea ice model‐data synthesis, is presented. ASTE_R1 has a nominal resolution of 1/3° and spans the period 2002–2017. The fit of the model to an extensive (O(109)) set of satellite and in situ observations was achieved through adjoint‐based nonlinear least squares optimization. The improvement of the solution compared to an unconstrained simulation is reflected in misfit reductions of 77% for Argo, 50% for satellite sea surface height, 58% for the Fram Strait mooring, 65% for Ice Tethered Profilers, and 83% for sea ice extent. Exact dynamical and kinematic consistency is a key advantage of ASTE_R1, distinguishing the state estimate from existing ocean reanalyses. Through strict adherence to conservation laws, all sources and sinks within ASTE_R1 can be accounted for, permitting meaningful analysis of closed budgets at the grid‐scale, such as contributions of horizontal and vertical convergence to the tendencies of heat and salt. ASTE_R1 thus serves as the biggest effort undertaken to date of producing a specialized Arctic ocean‐ice estimate over the 21st century. Transports of volume, heat, and freshwater are consistent with published observation‐based estimates across important Arctic Mediterranean gateways. Interannual variability and low frequency trends of freshwater and heat content are well represented in the Barents Sea, western Arctic halocline, and east subpolar North Atlantic. Systematic biases remain in ASTE_R1, including a warm bias in the Atlantic Water layer in the Arctic and deficient freshwater inputs from rivers and Greenland discharge.

Some of the recent changes in the observed Arctic Ocean heat content have been linked to pulsed warming of the Atlantic Water (AW) inflow (Muilwijk et al., 2018;Polyakov et al., 2017) and can be traced back upstream into the subpolar North Atlantic (SPNA; e.g., Årthun & Eldevik, 2016). Given the importance of Arctic changes and their interaction with the SPNA to the global climate system (Carmack et al., 2016), investigations of mechanisms setting the time-mean and evolving state of the Arctic Ocean and exchanges with surrounding ocean basins must be supported by basin-scale estimates of the ocean-sea ice state.
Historically, due to extremely sparse observations, efforts to construct decadal Arctic-focused gridded datasets have been hampered. Realistic simulation of the Arctic and sub-Arctic ocean-sea ice state has remained difficult due to highly uncertain initial conditions. Beginning in the early 2000s, increased availability of in situ observations of sub-surface ocean hydrography and of oceanic transports across Arctic gateways has improved our understanding of key processes, including interior eddy activity and mixing (Bebieva & Timmermans, 2016;Cole et al., 2014;Timmermans et al., 2012;Zhao et al., 2016), and the transformation and redistribution of watermasses (Proshutinsky et al., 2009;Pnyushkov et al., 2015;Polyakov et al., 2017;Rabe et al., 2014;Timmermans & Jayne, 2016;Timmermans et al., 2018;von Appen, Hatterman et al., 2015). Over the same time period, new satellite altimetry (Kwok & Morison, 2016), gravimetry (Peralta-Ferriz et al., 2014), and sea ice observations have allowed a more accurate estimate of Ekman transport (Meneghello et al., 2018) and inventory of freshwater in the Western Arctic (Proshutinsky et al., , 2020. In parallel with increased observational coverage, great progress has also been made using theoretical and modeling frameworks to advance our understanding of Arctic Ocean dynamics, for example, elucidating the importance of eddies in gyre equilibration (Manucharyan & Isachsen, 2019;Meneghello et al., 2017) and vertical heat redistribution (Polyakov et al., 2017). Despite this progress, confident assessment of the time-mean state, interannual variability and identification of robust decadal trends remains challenging (Balmaseda et al., 2015;Timmermans & Marshall, 2020) due to multiple factors (Docquier et al., 2019;Holloway et al., 2007;Ilicak et al., 2016;Q. Wang et al., 2016aQ. Wang et al., , 2016b. The most important amongst these factors is the lack of direct observations throughout the full water column, including at the air-ice-ocean interface, in and just below the mixed layer, along the AW boundary current pathway, and at the shelf-basin regions that connect the dynamics of this energetic current and the relatively quiescent Arctic Ocean interior (Timmermans & Marshall, 2020).
To fill these gaps, the community has constructed climatologies (e.g., WOA13 version 2 and WOA18, Locarnini et al., 2018;Zweng et al., 2018) and data-model syntheses (Carton et al., 2019;Stammer et al., 2016;Uotila et al., 2019) which are assumed to have higher fidelity as the repository of incorporated data grows. The improved fit between the latest climatology and existing observations is far superior to that seen in older climatologies. For example, in the Western Arctic interior, Ice Tethered Profilers (ITP) consistently report warmer temperatures ( Figure 1) than provided by both the Polar Hydrographic Climatology (PHC, Steele et al., 2001) and the World Ocean Atlas 2009 (WOA09, Antonov et al., 2010;Locarnini et al., 2010), but are in close agreement with WOA18. Here, it is important to note that this close agreement at the time/ location of data acquisition is built into the majority of these climatologies and other existing Arctic model-data syntheses. These products are constructed using statistical methods such as optimal interpolation (e.g., PHC, WOA), 3D-Var, or sequential 4D-Var with short assimilation windows (Carton et al., 2019;Mu et al., 2018;Stammer et al., 2016;Uotila et al., 2019). The advantage of these methods is that the synthesis ensures a local fit to available observations ( Figure 1, Carton et al., 2019). Away from observed locations, however, the interpolator relies on incomplete, unavailable or unobtainable information. Missing values are, for example, determined via spatial/temporal correlations, potentially derived from regions/times of very different dynamics. By construction, high frequency variability cannot be fully accounted for and as a result spectral agreement with observations can be poor (Verdy et al., 2017). Importantly, this type of interpolation-and that used in 3D-Var or sequential 4D-Var-can introduce artificial sources/sinks (e.g., of mass, enthalpy, and momentum, Griffies et al., 2014;Stammer et al., 2016;Wunsch & Heimbach, 2013), which make a large contribution to the total energy budget (Balmaseda et al., 2015). This violation of basic NGUYEN ET AL. 10.1029/2020MS002398 2 of 49 conservation principles has been shown to obfuscate the use of these products for robust identification and attribution of change, creating spurious trends (Bengtsson et al., 2004), and triggering artificial loss of balance (Pilo et al., 2018), resulting in adjustments that may propagate and amplify to corrupt the large scale solution (Sivareddy et al., 2017).
To lend additional support to studies of the Arctic ocean-sea ice system over the early 21st century we have developed a new model-data synthesis utilizing the non-linear inverse modeling framework developed within the consortium for Estimating the Circulation and Climate of the Ocean (ECCO, Heimbach et al., 2019;Stammer et al., 2002;Wunsch & Heimbach, 2007). The use of the primitive equations as a dynamical interpolator distinguishes our effort from purely statistical approaches. The inversion consists of an iterative, gradient-based minimization of a least squares model-data misfit function. Unlike most reanalysis products that are based on sequential data assimilation, only independent, uncertain input variables, that is, initial conditions, surface boundary conditions and model parameters are adjusted. No periodic analysis increments during the estimation period that would incur artificial sources or sinks are permitted. Through strict adherence to conservation laws, all sources and sinks within the state estimate can be accounted for over the full estimation period, permitting meaningful analysis of closed budgets (Buckley et al., 2014;Piecuch & Ponte, 2012).
Our work builds upon extensive prior efforts of the ECCO community to produce optimal (in a least squares sense) kinematically and dynamically consistent data-constrained estimates of the ocean state across the globe and in various regional domains. Among the publicly available ECCO state estimates is ECCO Version 4 Release 3 (ECCOv4r3, Forget, Campin et al., 2015;, which has been constrained to satellite and in situ data (including Argo and elephant seal data) outside of the Arctic, ITP data in the Arctic, and other mooring data at important Arctic gateways. The ECCOv4 releases have been widely used, with applications including investigation of global vertical heat and salt redistribution (Liang et al., 2017;Liu et al., 2019), heat budgets in the North Atlantic (Buckley et al., 2014(Buckley et al., , 2015Foukal & Lozier, 2018;Piecuch et al., 2017) and the Nordic Seas (Asbjørnsen et al., 2019), high-latitude freshwater budgets (Tesdal & Haine, 2020), and sea level change (Piecuch & Ponte, 2013).
The state-estimation procedure entails reducing the total time-and space-integrated model-data misfit. Since ECCOv4r3 is a global solution, reduction of the relatively well-sampled misfit at lower latitudes dominates the production of this solution. As a result, ECCOv4r3 possesses notable biases in the Arctic (Carton NGUYEN ET AL.   Tesdal & Haine, 2020), including a strong anticyclonic circumpolar circulation of Atlantic Water ( Figure 2). Furthermore, the ECCOv4r3 horizontal grid spacing of 40-45 km is well above the Rossby deformation radius in the Arctic and Nordic Seas (Nurser & Bacon, 2014). This has motivated a dedicated effort to build a higher resolution regional state estimate for use in Arctic inter-annual to decadal climate research, covering the early 21st century, culminating in the Arctic Subpolar gyre sTate Estimate (ASTE).
Here we describe the first release of ASTE (ASTE_R1), providing an estimate of the ocean-sea ice state for the period 2002-2017. We describe the model configuration, observational constraints and the state estimation machinery (Section 2) and present the model-data misfit reduction (Section 3). We then compare our estimates of volume, heat and freshwater transports through important Arctic gateways with those in the existing literature as well as present an analysis of ASTE_R1 heat and freshwater budgets for the Arctic Ocean, Greenland-Iceland-Norwegian (GIN) Seas and subpolar North Atlantic (section 4). In Section 5, we examine how an improved fit is achieved, identifying key adjustments of our independent control variables, and review remaining issues in ASTE_R1. In Section 6 we summarize key findings and discuss future directions.

Model Description
The coupled ocean-sea ice model underlying the estimation framework is an evolved version of the Massachusetts Institute of Technology general circulation model (MITgcm; Adcroft et al., 2018;Marshall et al., 1997). The model solves the primitive equations in rescaled z* coordinates  with a full non-linear free surface . The dynamic-thermodynamic sea ice model is an evolved version of Menemenlis et al. (2005); Losch et al. (2010) and Heimbach et al. (2010). Eddy-induced tracer mixing and transports along isopycnal surfaces are parameterized following Redi (1982) and Gent and McWilliams (1990).
The model uses a finite-volume discretization in a so-called "latitude-longitude-polar-cap" grid configuration (LLC grid, Forget, Campin et al., 2015). The LLC grid is topologically equivalent to a cubed-sphere grid , but reverts to a regular latitude-longitude grid equatorward of ∼57°N. The computational cost associated with solving the non-linear optimization problem for eddy-resolving simulations, which would require resolutions well below 4-15 km for the Arctic Mediterranean (Nurser & Bacon, 2014), is prohibitively high. As a compromise, ASTE is based on the medium-resolution LLC-270 grid, providing NGUYEN ET AL.  a nominal grid spacing of 1/3°, which corresponds to ∼22 km in the North Atlantic, ∼16 km in the Nordic Seas, and ∼14 km in the high Arctic interior (Figure 3).
The ASTE domain covers the entire Atlantic northward of 32.5°S, the entire Arctic and its surrounding seas (Labrador, Nordic, Barents, Bering north of 47.5°N) and the Canadian Archipelago. The model has 50 unevenly spaced vertical height levels; thicknesses range from 10 m at the surface to 500 m at 5,000 m depth. The 10 m thickness at the surface cannot fully resolve surface boundary layer processes or the shallowest summer mixed layer of ∼5 m, but is deemed sufficient for capturing the 10-100 m seasonal MLD in the Arctic (Bigdeli et al., 2017;Peralta-Ferriz & Woodgate, 2015;Rudels, 2015;Rudels et al., 2004) and is a reasonable choice given the size and expense of our computations. Partial cells  are used to improve the representation of topography. The domain has boundaries at 35°S in the South Atlantic, 48.6°N in the Pacific, and at the Gibraltar Strait. Rationales for choosing a full Atlantic-Arctic domain for ASTE-rather than limiting it to the Arctic Mediterranean-are to extend the applicability of the solution to investigation of latitudinal connectivity between Atlantic and Arctic variability on decadal timescales, and to displace the imposed open boundary conditions far from the region of key interest.
We prescribe lateral open boundary conditions from the global ECCOv4r3 solution, which has been shown to be in good agreement with large-scale constraints from satellite and in situ data (including Argo). The bathymetry is a merged version of Smith and Sandwell (1997), version 14.1, below 60°N and the international bathymetric chart of the Arctic Ocean (Jakobsson et al., 2012) above 60°N, blended over a range of ±100 km about this latitude. Special attention was paid to remove abrupt jumps over the merged region. Model depths within important canyons (e.g., Barrow) and across important gateways (e.g., Florida Straits, Greenland-Iceland-Faroe-Scotland ridge, Aleutian Islands chain, Gibraltar Strait) were enforced to be consistent with observations in order to realistically simulate key transports and regional circulations.
Atmospheric forcing is applied via bulk formulae (Large & Yeager, 2008) over the open ocean, with the initial estimate of the atmospheric state variables from JRA-55 (Kobayashi et al., 2015). We considered taking ERA-Interim (Dee et al., 2011)-employed by ECCOv4r3 (Forget, Campin et al., 2015;)-as our first guess. However, this product has a well documented warm bias of up to 2°C in the Arctic (Beesley et al., 2000;Freville et al., 2014;Jakobson et al., 2012;Lupkes et al., 2010)   Lawrence, and all channels in the Canadian Arctic Archipelago except Nares and Barrow Straits, are masked. Depths of several important channels, including the Barrow Canyon, Greenland-Iceland-Faroe-Scotland Ridge and the Florida Straits, were carefully inspected to ensure transports consistent with published observations. the Japanese Reanalysis (JRA-25) without the need to increase sea ice albedos above their observed values. For this reason, the updated 3-hourly, higher-resolution JRA-55 was chosen as the initial surface boundary forcing. Monthly mean estuarine fluxes of freshwater are based on the Regional, Electronic, Hydrographic Data Network for the Arctic Region (R-ArcticNET) data set (Lammers & Shiklomanov, 2001;Shiklomanov et al., 2006). Figure 2 the ECCOv4r3 solution does not exhibit the cyclonic circulation of Atlantic water in the Arctic that is inferred from hydrographic observations (Rudels, 2012). For this reason, we elected to initialize from alternative products. Table 1 summarizes our first-guess model input parameters for sea ice, ocean mixing and momentum dissipation, along with our choice of ocean-sea ice state to initialize the unconstrained simulation. This run serves as iteration 0 of the optimization and will be referred to as it0 for the remainder of the study. Our selection is informed by existing observation/model-based estimates. Importantly, sea ice albedos and drag coefficients are chosen within the range of observed and previously optimized estimates (Nguyen et al. 2011).

As shown in
The three-dimensional parametric horizontal stirring fields for temperature and salinity are based on typical values used in the literature (Pradal and Gnanadesikan (2014); Campin, 2014 (pers. comm.)) with consideration for where the ASTE grid resolves the baroclinic deformation radius as follows. The vertical background diffusivity z  was set based on typical values at latitudes below 79°N of ∼10 −5 m 2 /s and limited observed and modeled ranges of 10 −7 -10 −6 m 2 /s at high latitudes (Cole et al., 2014;Fer, 2014;Nguyen et al., 2011;Padman & Dillon, 1988;Sirevaag & Fer, 2012;Zhang & Steele, 2007). Vertical diffusivities are enhanced by a factor of 10 near the sea floor to mimic lee wave-driven mixing (Mashayek et al., 2017;Toole, 2007). Horizontal dissipation is applied as a combination of biharmonic Leith and Laplacian viscosity (Fox-Kemper & Menemenlis, 2008;Griffies, 2004). At lower latitudes, where eddy effects are better resolved, we follow the formulation of Leith (1996) to represent the direct enstrophy cascade at mesoscales. Within the attached Gulf Stream, a higher Laplacian viscosity is initially required to reduce the Reynolds number and prevent premature separation (Chassignet & Garraffo, 2001;Chassignet & Marshall, 2013;Dengg, 1993). In the Arctic Mediterranean, where the deformation radius is 4-10 km, an ad-hoc combination of biharmonic Leith and Laplacian viscosity is used to ensure consistency of inflow velocity at Fram Strait and an approximate cyclonic circumpolar AW circulation inside the Arctic (Jochum et al., 2008, see Table 1). The model is spun up for 6 years using repeated year 2002 atmospheric forcing and open boundary conditions (Table 1). The ocean, sea ice and snow states at the end of this 6 years spin up became the initial condition for the unconstrained it0 in the optimization procedure described next.

State Estimation Framework
ASTE is formally fit to observations through a gradient-based iterative least-square minimization of the model-data misfit function that takes into account data and model parameter uncertainties . The gradient with respect to a high-dimensional space of uncertain input variables, the "controls," is obtained via the adjoint of the model, derived by means of algorithmic differentiation (AD; Giering et al., 2005;Heimbach et al., 2005). The model-data misfit (or "cost") function is defined as (Wunsch & Heimbach, 2007): where time t ∈ [t 0 , t f ], t 0 and t f are the initial and final time, and Δt the time-stepping of the forward model. y(t) is the observation vector and x(t) the state vector containing the model ocean (e.g., temperature, salinity, velocities, sea surface height) and sea ice variables (e.g., concentration, ice and snow thickness, velocities) at all grid points (Wunsch & Heimbach, 2007). The combined initial model state x 0 and input parameter adjustments u(t) collectively comprise the control vector Ω ∋{x 0 , u(t)}. E is the operator mapping NGUYEN ET AL.
10.1029/2020MS002398 6 of 49 the state variables to the observations. The model-data misfit y(t)−E(t)x(t) is weighted by the inverse error covariance matrix R(t). This accounts for both observational uncertainty and model representation error, where the latter considers the extent to which real variability cannot be represented at the chosen model resolution . B(t 0 ) and Q(t) are error covariances of x 0 and u(t), respectively. Full knowledge of R, B, and Q is often unattainable (Wunsch & Heimbach, 2007). As a result, the misfit y(t)−E(t)x(t) and variables δx 0 = x 0 − x(t 0 ) and u(t) are often assumed Gaussian, with zero means and standard deviations whose squares fill the diagonal entries of their respective covariance matrices (Wunsch & Heimbach, 2007). In the absence of better information, we resort to the simplified representation of the error covariances, consistent with existing state estimation efforts (e.g., Forget, Campin et al., 2015;Mazloff, Heimbach, & Wunsch, 2010). These error estimates play an important role in any least squares optimization (both ECCO-related and other data assimilation efforts), and their improved estimation is itself an important area of ongoing research (Wunsch, 2018). We will discuss these further below.
There are three distinct contributions to the misfit cost function, Equation 1. The first term describes the normalized model-data squared misfit to be minimized. This term sums weighted contributions from all observational data considered. The second term penalizes deviation of the initial state x(t 0 ) from the initial guess x 0 (Table 1). Similarly, the third term describes moderation of input parameter adjustments u(t) so that the adjustment amplitude does not far exceed the uncertainties. The adjoint (or Lagrange multiplier) method consists of augmenting the cost function (Equation 1) to a Lagrangian function  by adding an additional term that enforces the strict adherence of the solution to the model equations. In this manner, the constrained optimization problem (find extrema of J subject the constraint that the model equations be fulfilled exactly) is converted into an unconstrained problem of finding stationary points of the Lagrangian (Wunsch & Heimbach, 2007).
The optimization problem is solved via gradient-based optimization, in which the gradient of the cost function with respect to the control variables informs an iterative minimization algorithm. In our case, this is the quasi-Newton method following Gilbert and Lemaréchal (1989). Once the cost function J is defined, beneficial control adjustments that reduce the misfit are informed by the gradient ∇ Ω J. This gradient can be efficiently computed for very high-dimensional control spaces using the adjoint model (Wunsch & Heimbach, 2007). These adjusted controls are then used in a new integration of the forward model for the full period (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017), during which model-data misfits are recomputed. At the end of this forward integration, contributions to the costfunction are accumulated, the adjoint model is integrated, and the gradient information is re-computed informing updated control adjustments for the next integration of the forward model. The optimization thus proceeds in an iterative manner, whereby each iteration entails execution of both the forward and adjoint model, providing updated control adjustments to obtain further reduction of the total model-data misfit in successive iterations. The optimization is continued until little further misfit reduction is achieved between successive iterations. This is expected when the state estimate is in agreement with the observations within the error R(t) (expressed, e.g., in terms of a χ 2 distribution of the squared normalized misfit residuals).
The space of control variables for ASTE, {x 0 , u}∈Ω, comprises the 3D hydrographic initial conditions, potential temperature and salinity (θ 0 , S 0 ), the time-varying 2D surface atmospheric state variables, spatially varying but temporally invariant model coefficients of vertical diffusivity ( z  ) and parameterized eddy activity (  , gm   ), denoting the strength of eddy-induced isopycnal diffusivity and potential energy transfer, respectively (Forget, Ferreria et al., 2015). The atmospheric state control variables are 2 m air temperature, T air , specific humidity, q air , downward short-and long-wave radiation, R sw , R lw , precipitation, P, and 10 m winds u w , v w . Although runoff and evaporation are not control variables, in practice they project onto the precipitation sensitivities, interpreted as the linear combination of net surface freshwater fluxes (evaporation minus precipitation minus runoff, E-P-R).
To ensure the adjustments are physically reasonable, a-priori uncertainties (i.e., the square-roots of the diagonal terms of B and Q) are estimated following Forget and Wunsch (2007); Fenty and Heimbach (2013a); Fukumori, Fenty et al. (2018) for oceanic hydrography and Chaudhuri et al. (2013Chaudhuri et al. ( , 2014 for atmospheric forcing. Whilst the estimate from Forget and Wunsch (2007) and Fenty and Heimbach (2013a) Chaudhuri et al. (2013Chaudhuri et al. ( , 2014 are based on the spread between atmospheric reanalysis products, which is particularly large over the Arctic.  The vector y(t) contains as many available ocean and sea ice observations as we were able to access. The observational backbone of ASTE_R1 includes the standard ECCOv4r3 suite (Table A1) of in situ and remotely sensed ocean data: temperature and salinity profiles from Argo, GO-SHIP and other research cruises, instrumented pinnipeds, gliders, and moorings, and ice-tethered profilers; ocean bottom pressure anomalies from GRACE (Watkins et al., 2015;Wiese et al., 2018); sea surface height from Ocean Surface Topography Mission/Jason 2 and Jason 3 (Zlotnicki et al., 2019); Mean Dynamic Topography DTU13 (Andersen et al., 2015); and infrared and microwave-derived sea surface temperature (JPL_MUR_MEaSUREs_Project, 2015). For details on how the data and their uncertainties were obtained and prepared we refer the readers to Fukumori, Fenty et al. (2018). In addition to the ECCOv4r3 suite, the data are augmented by updated high latitude in situ profiles, ship-based CTD, and mooring observations at important Arctic gateways and in the Arctic interior (see Table 2, Figure 4).
The estimation period chosen for ASTE, 2002-2017, leverages the increase in satellite (GRACE, ICESat-1/2, CryoSat-2) and in situ (ITP) observations in the Arctic, as well as the beginning of the quasi-global Argo float deployment. In total, ∼1.2 × 10 9 observations were employed to constrain distinct aspects of the modeled ocean and sea-ice state, culminating in the optimized ASTE_R1 solution. Key among these are satellite-based observations of sea level anomalies (SLA) to aid removal of the precipitation bias in JRA-55, Argo and lower latitude CTD to improve surface and sub-surface hydrography in the North Atlantic and Nordic Seas, and a suite of moorings in the Arctic. This suite includes the Fram Strait mooring array to constrain the boundary current strength and heat flux from the Nordic-Seas into the Arctic, and the combined ITP and Beaufort Gyre moorings to constrain the Canada Basin hydrography. Finally, OSISaf daily sea ice concentration was essential for constraining the ice edge and upper ocean hydrography in the Arctic and its surrounding marginal seas.
The error R associated with the observations y(t) is the combined data uncertainty and model representation error. For hydrographic data, the derivation is as described above for the a-priori uncertainties B. For satellite data, errors are a combination of the corresponding satellite mission's provided uncertainty and model representation errors, as described in Fukumori, Fenty et al. (2018). These representation errors were derived from the data variance within (and weighted by the area of) ASTE horizontal grid box. They generally exceed the stated mission uncertainty. As seen in Equation 1, R plays an important role in weighting the individual model-data misfit terms. Careful assessment of R is thus required to ensure an appropriate and balanced contribution of the diverse datasets to the total J.
The practical implementation of Equation 1 follows that described in Forget, Campin et al. (2015). Several approximations to parameterization in the adjoint model were made to ensure stable behavior. Maximum isopycnal slopes are limited in the GM/Redi parameterization, the vertical mixing scheme (K-Profile Parameterization, Large et al., 1994) is omitted in the adjoint, and increased horizontal and vertical momentum dissipation are employed in the adjoint to suppress fast growth of unstable sensitivity.
The full sea ice adjoint, as described in Fenty and Heimbach (2013a); Fenty et al. (2015), was not used in this study (nor in ECCOv4), due to persistent instability issues. In its place, the sea ice concentration model-data misfit is used to relate air-sea fluxes to the enthalpy of the integrated surface ocean-sea ice system as follows.
Where the model has an excess/deficiency of sea ice, extra heat is added to, or removed from the system to bring the sea surface to above or below the freezing temperature. In these two cases, the pseudo-sea ice cost function contributions J seaice_conc_ [ex,de] are in enthalpy units rather than normalized model-data misfits. Normalization is chosen to obtain amplitudes comparable to other model-data misfits J i contributing to the total cost function J, so that these terms play an active role in the optimization. Lastly, convergence-if achievable for these two pseudo-sea ice costs-is when they approach zero and not unity.
After 62 iterations, a substantial reduction in model-data misfit has been achieved compared to the unconstrained simulation, such that the solution is deemed suitable as ASTE first release. The initial conditions of the optimized state, ASTE Release 1 (ASTE_R1), are derived by adding the adjustments [Δθ, ΔS] i62 to the first guess fields [θ, S] i0 . The same holds for the optimized mixing fields and surface atmospheric state (see Table 1). Adjustments to the uncertain control variables obtained as a result of the gradient-based optimization enable the improvement in the model fit to observations while retaining dynamical consistency. Figure 5 shows the uncertainty and adjustments for four of the seven surface atmospheric state variables, NGUYEN ET AL.   Chaudhuri et al. (2014), shows some of the largest disagreements amongst the atmospheric reanalyses to be in the Arctic (Figure 5a1-5a4). The percentile (pctl) thresholds indicate that for these four fields the adjustments are within the uncertainty. Overall, the 99-pctl adjustments are within 2σ for all time-dependent atmospheric variables except downward shortwave where it is within 3σ.
The full monthly mean state of ASTE_R1 is distributed via the ECCO & ASTE data portal at the Texas Advanced Computing Center (TACC). In addition, the model configuration, required input fields and code are distributed to enable reruns (see Appendix A). Since the focus of ASTE is on the North Atlantic and the Arctic Ocean, we restrict both the discussion presented below and the distributed ASTE_R1 fields to latitudes above 10°N. As for ECCOv4r3, the mass, salt, and heat budgets in ASTE_R1 are accurately closed when computed using the distributed standard ECCO diagnostics that we provide. In Appendix B, we show how lateral transports may be accurately computed and provide estimates for the errors incurred in offline calculations using the ASTE_R1 monthly mean diagnostics.

Making Meaningful Model-Observation Comparisons
A meaningful assessment of ASTE_R1 through comparison with observations is nontrivial and requires careful consideration. One of the biggest challenges is properly accounting for the sparse spatio-temporal sampling and the potential for aliasing. For example, measurements might only be taken at a discrete location (e.g., a mooring) or along-track (i.e., with high along-track coverage and drastically lower resolution in the cross-track direction) or only during summer months (e.g., ship-based CTD). "Averages" of these measurements (e.g., average Argo or ITP data over 1 month or 1 year) incur aliasing in both space and time as well as potential spatial or seasonal biases. "Averages" of ASTE_R1 outputs at the smallest spatial scale (grid cell size), on the other hand, are over a spatial area of ∼200 km 2 . Unless observations are well sampled over this grid-area, a direct comparison between observations and ASTE_R1 can be problematic. Furthermore, ASTE_R1 does not resolve eddies in the Arctic and GIN Seas. As a result, we should neither expect nor demand a perfect fit to discrete (in space/time) measurements. As is common in data assimilation (Janjić et al., 2017) the ECCO framework utilizes "representation errors"-described above-in the weighting of the model-data misfits, Equation 1, to safeguard against over-fitting and facilitate more meaningful model-data comparison (Wunsch & Heimbach, 2007). However, these representation errors are themselves highly uncertain, often relying on unconstrained high-resolution model runs from which they are inferred (see Nguyen, Heimbach et al., 2020, for a more detailed discussion).
NGUYEN ET AL.
10.1029/2020MS002398 11 of 49 In Section 3 we present both normalized misfit reductions as well as comparisons of dimensional transports and heat/freshwater contents. For the dimensional quantities, we encounter several potential challenges related to resolution and bias issues, which we briefly discuss in the following. A serious challenge stems from the need to compare watermasses in the presence of hydrographic biases. From observations, watermasses are often defined in temperature, salinity, and density (T, S, and σ) space with tight thresholds/bounds reflecting the measurement precision (e.g., to the first or second decimal place). These definitions can be problematic to adopt in ASTE_R1, where we are averaging over grid-cell areas of ∼200 km 2 and thicknesses of 10-500 m. Furthermore, in some regions the model representation errors may be up to an order of magnitude larger than measurement precision. In these regions, the normalized misfit can be within acceptable range, but ASTE_R1 can still possess notable absolute (T, S, and σ) biases if the representation errors are large. For this reason, a watermass is likely to exist in ASTE_R1 but with modified thresholds/bounds. Where appropriate, we analyzed ASTE_R1 carefully in T, S, and σ space to identify suitable classifications for calculation of watermass transports. Details on the modified bounds are provided in Appendix C.
A second challenge is related to the region over which derived quantities are computed. In cases where these regions are defined with geographic bounds based on availability of observations rather than dynamical regimes, the equivalent derived quantities in ASTE_R1 can be highly sensitive to small shifts in bounding region, especially when the grid resolution and uncertainties in the control input parameters (e.g., forcing, internal mixing) are taken into account. For this reason, we also explore the sensitivity of area/volume integrals to choice of geographical bounds in Appendix C.
Lastly, comparison of (dimensional) integrated transports can be problematic due to spatial sampling issues and representation error, preventing precise estimation of narrow boundary currents in ASTE_R1. An example is at Fram Strait, where the ASTE_R1 grid cannot resolve the e-folding scale of the West Spitsbergen NGUYEN ET AL.  . In situ observations used to constrain ASTE. Red squares with "x" are additional OSNAP mooring data, used for independent evaluation but not part of the cost function. Figure 5. The atmospheric forcing field (a1-a4) uncertainty σ and (b1-b4) the 50-percentile and (c1-c4) 90-percentile normalized adjustment magnitudes |δ|/σ. The uncertainty fields have units given above the color scale for a1-a4. Note that the reciprocal of the squared value of these uncertainties are entries in the weight matrix Q. The normalized adjustment magnitudes are dimensionless.
Current (Beszczynska-Möller et al., 2012). Enforcing fit to the observed mooring velocity would likely result in an overestimation of the net inflow volume transport here. In ASTE_R1, velocities at gateways were not employed as active constraints but were used for offline assessment of the derived transports. Ultimately, however, the spacing between discrete moorings offers incomplete information on the total volume transports across a given gateway, and existing observation-based estimates generally require various assumptions on spatial/temporal correlations in order to interpolate between the mooring measurements. As a result, our direct comparisons of ASTE_R1 and observation-based transports presented below seeks consistency in terms of sign and order of magnitude rather than exact agreement of amplitude. This is especially true for assessment of ASTE_R1 heat and freshwater transports, computed relative to the wide range of reference values used in the literature.

Model-Data Misfit Reduction and Residuals
In what follows, we will assess the ASTE_R1 solution in the context of existing observation-based estimates of the circulation and hydrography in the Arctic. We first compare it0 and ASTE_R1 using the online and offline cost metrics described in Section 2 and listed in Table 3, and summarize the reduction in the integrated model-data misfits and costs achieved in the production of ASTE_R1. We then expand this discussion, considering the ASTE_R1 fit to constraints in the Arctic, GIN Seas, and Subpolar North Atlantic (Sections 3.1-3.3). Note that assimilation aids-but by no means guarantees-model-data consistency due to errors and/or deficiencies in the data, model, and/or state estimation framework. This point will be revisited in our discussion in Section 5. We refer to "misfit" as the dimensional model minus data difference, "normalized misfit" as misfit scaled by the respective uncertainty (dimensionless), and "normalized cost" or "cost" as the square of the normalized misfit (dimensionless). The overall cost reductions in ASTE_R1 have been grouped into several categories as shown in Figure 6 and summarized in Table 3.

Sea Ice
Improved representation of sea ice extent in ASTE_R1 (compared to the unconstrained simulation) is indicated by a significant reduction of by 53% and 42%, respectively (Table 3, Figure 6). Fenty and Heimbach (2013b) showed that these improvements can be effectively achieved through small adjustments of atmospheric controls, within their uncertainty range. These improvements are independently confirmed by the reduction in offline misfits for sea ice area (J seaice_area15 , 64%) and extent (J seaice_extent15 , 83%) ( Table 3). The largest improvements occur in the seasonal ice zones e.g., Greenland and Barents Seas and Southern Beaufort Gyre (Figure 7a) associated with a systematic decrease in total simulated area/extent without alteration of the seasonal cycle ( Figure 7b).
The improved sea ice edge representation in ASTE_R1 (Figures 7e and 7h) is accompanied by a reduction in the offline misfits for sea ice velocities (J seaice_vel , Table 3), primarily in locations where nonzero ice velocities in it0 were accompanied by observations of zero ice concentration and vice versa. Unlike velocity, however, the sea ice thickness costs J seaice_thickness did not decrease ( Note. The quantities listed above the triple horizontal lines contribute directly to the total J in Equation 1, i.e., contribute to the gradient-based minimization, whereas those listed below the thick horizontal lines (J (o) ) are purely diagnostic, i.e., are used only for offline assessment and do not influence the optimization. The offline sea ice area and extent (both defined using the common 15% cutoff threshold) costs are normalized by the Arctic Mediterranean's area.

Table 3 Active and Offline Costs and Reductions in ASTE_R1 Compared to it0
adjoint does not contain physics relating ice thickness to the atmospheric forcing or ocean interaction from below. We will return to this in section 6.

Fram Strait
The dynamics in the vicinity of Fram Strait are highly complex, governed by strong air-ice-ocean interaction, vigorous generation of eddies associated with highly sheared boundary currents and their recir- In view of the importance of AW for the wider Arctic region, we paid particular attention to skillfully represent AW inflow at Fram Strait as follows. Early in the development of ASTE_R1 and prior to the gradient-based optimization, we compared the simulated volume transport to daily average moored velocity (available for the years 2002-2011) at various depths along the entire array , 2012. The model viscosity was prescribed (Section 2.1, Table 1) to ensure we obtained a representative volume transport across the strait. During the iterative optimization, the mooring temperature and salinity were included in the active costfunction (J) to directly constrain T/S at the strait. Via these steps, indirect constraint of the tracer transports (i.e., v * T and v * S) and their constituents (e.g., "inflow/outflow of AW") at Fram Strait was achieved, as shown by comparison to additional data for the years 2012-2017 (von Appen,  that were withheld from the optimization for offline evaluation. Figure 8 shows the unconstrained it0 and ASTE_R1 misfits to moored T, S, and northward velocity, as a function of longitude. For the region occupied by AW inflow, the normalized misfit in temperature is reduced by 71% compared to the unconstrained it0 simulation (Figure 8a). On the continental shelf, west of 4°W, ASTE_R1 has a warm bias compared to the observations, which resulted in higher misfits here (red bar between longitude 8.1°W and 4.1°W in Figure 8a). The reduction in misfits for inflow of the AW, however, is more important for the large-scale Arctic hydrography, as AW passing through this important gateway propagates along the entire boundary of the Eastern Arctic and into the Canada Basin. Here its properties can be compared to ITP data, which serve as the main constraint on subsurface T/S over the entire pathway from the Fram Strait (see Section 3.2.3). The misfit reduction can be seen for one example mooring at approximately 8°E (Figures 8d-8f) at multiple depths. The significant improvement in salinity at the surface (dashed blue in Figure 8e) is related to the improved ice edge (see also Figures 7c-7e). Another significant improvement is in the AW core temperature at depth ∼250 m (solid blue and red lines in Figure 8d).
NGUYEN ET AL.

10.1029/2020MS002398
14 of 49 Figure 6. Aggregated cost reductions calculated from key data sets that were used in the optimization. The numbers listed above each data set are the percentage of cost reduction in ASTE_R1 compared to it0. The magnitudes of the two pseudo-sea ice costs are indicated on the right abscissa.

Canada Basin Hydrography
Once the AW, via the West Spitsbergen Current, crosses Fram Strait and traverses the Eastern Arctic along the continental slope Polyakov et al., 2017;Rudels, 2015), the watermass properties (e.g., current strength and direction, density, temperature) are not as well constrained due to extreme data paucity in the Eastern Arctic. In particular, along the boundary current path (shoreward of the red contour in Figure 9a), only 1% of the total 2004-2016 ITP data are acquired within the Nansen Basin; only 4.5% are acquired in the combined Amundsen and Makarov Basins (Figures 9b-9d). The majority of the ITP data (71%) are within the Canada Basin ( Figure 9e). Figure 10 shows the misfit reduction between the unconstrained it0 and ASTE_R1 as a function of basins and depths. The reduction is throughout the upper 800 m of the water column. Seawater density in the Arctic is primarily controlled by salinity, whereas temperature behaves more like a passive tracer and can be more flexibly impacted by the optimization procedure. As a result, the reduction in ASTE_R1 temperature misfits greatly exceeds the reduction in salinity misfits in the Arctic, with the most notable improvement occurring within the AW core (120-450 m, Figures 10a-10c). The remaining notable temperature misfits in ASTE_R1 are at depths occupied by the mixed layer (10-65 m) and below the AW core (450-800 m, Figure 10c). Large salinity misfits also persist in the mixed layer and in the halocline (120-250 m, Figures 10d-10f).
NGUYEN ET AL.

10.1029/2020MS002398
15 of 49 The optimization has, nevertheless, significantly narrowed the misfit distribution, eliminating the largest amplitude biases throughout the water column in temperature and especially below 400 m depth in salinity. The overall reduction is 85% for temperature and 56% for salinity.
An example of how temperature misfits are reduced in the water column is shown in Figure 11 for ITP #55, whose trajectory began in the Canada Basin interior (red circle in Figure 11a) and ended at the slopes of the Chukchi Plateau (green square). In the observations several watermasses can be seen, including the surface cold layer above ∼30 m, warm Pacific Summer Water (PSW) at ∼40-100 m, Cold Halocline Waters at ∼110-250 m, and the Atlantic Water core at depths ∼300-750 m (Figure 11b). In the unconstrained it0, both the AW boundary current and the halocline are too warm, the AW layer is too thick, and the PSW is too cold. In ASTE_R1, closer consistency is obtained with temperature observations for all watermasses. We emphasize that we have not applied direct adjustments to the time-varying simulated ocean state to achieve this fit (i.e., no "analysis increments" were applied). Instead, it is achieved through adjustments of NGUYEN ET AL.
the control variables, that is, the initial hydrography in 2002, time-averaged internal mixing parameters, and surface atmospheric forcing. As a result, a "near-perfect" fit, such as that of the WOA18 hydrography to the mean ITP data seen in Figure   within the specified temperature and salinity uncertainties, with improved watermass representation for all ITP data (e.g., Figure 11), Beaufort Gyre Moorings, NABOS moorings, and Fram Strait moorings.

The Greenland Iceland Norwegian Seas
The Greenland-Iceland-Norwegian (GIN) Seas are defined here as bounded to the south by the Greenland-Scotland Ridge (GSR) and to the north and north east by the Fram Strait and the Barents Sea Opening, respectively (see Figure 12b). The sea ice near Fram Strait and along the East Greenland coast is seasonal and the largest misfits in it0 were due to excessive ice here, including the Odden ice tongue (Wadhams et al., 1996) reaching further to the east during winter months (Figures 7c, 7d and 12a). Surface winds and air temperature have been found to play an important role in controlling the eastern extent of the ice edge in this region (Germe et al., 2011;Moore et al., 2014). Adjustments of these atmospheric state variables during the optimization, within their specified uncertainties, drove a reduction in sea ice area (Figure 12a) to improve the model-data fit.
The Nordic Seas host the interaction of several important watermasses. Warm and salty Atlantic water enters across the GSR along three major branches, meeting locally modified water recirculating in the Lofoten, Greenland and Iceland Basins, and the southward flowing cold, fresh East Greenland Coastal Current (Hansen & Østerhus, 2000). This region is characterized by very weak stratification, resulting in a very small deformation radius of 4-7 km throughout the region (Nurser & Bacon, 2014), which further challenges realistic representation of watermass distribution in models (Drange et al., 2005;Heuzé & Årthun, 2019) and ASTE_R1. Nevertheless, the improvements obtained in ASTE_R1 are substantial, with overall reductions of ∼85% and 30% for salinity and temperature costs, respectively, through the 2,000 m water column (Figure 12c). The largest improvements are associated with reduction of a fresh bias in the upper 100 m (Figures 12b and 12c), across the Lofoten, Iceland, and Greenland basins, which are important regions for deep water formation.

The Subpolar Gyre and North Atlantic
Although the primary focus of the study is on the assessment of ASTE_ R1 in the Arctic Mediterranean, the North Atlantic Ocean serves as both the source of near surface heat and salt to the Arctic and the sink of dense deep water and surface freshwater from the Arctic Mediterranean, and so will be briefly assessed here.
One of the greatest challenges in modeling the North Atlantic is to correctly simulate the observed Gulf Stream pathway. Capturing a realistic Gulf Stream separation is non-trivial in z-level numerical models (Ezer, 2016;Chassignet & Xu, 2017    biharmonic and off-shore Leith viscosity as described in Section 2.1 was used to achieve an observationally consistent mean Florida Strait transport of ≈32 Sv (Baringer & Larsen, 2001;Johns et al., 2002) and a separation near Cape Hatteras. After separation, the Gulf Stream path can be approximately tracked using a proxy of the 15°C isotherm at 200 m depth (from the WOA13, Wolfe et al., 2019). Figure 13a shows this proxy of the Gulf Stream path for the years 2002-2017 in ASTE_R1 compared to that derived from WOA13.
The dynamical mechanisms underlying the transports of warm AW from the Gulf Stream extension to the subpolar North Atlantic (SPNA) and into the GIN Seas across the GSR is poorly understood and its representation in state-of-the-art models remains a great challenge (Heuzé & Årthun, 2019). Compared with Argo data, the eastern SPNA hydrography (south of the GSR) contains large biases in it0 (Figures 13b  and 13c) but is significantly improved in ASTE_R1, with a net reduction of misfit over the entire North Atlantic (north of 12°N) of ∼80% and ∼21% in salinity and temperature, respectively (Figure 13b). Closer inspection reveals that, just south of the GSR, in the Irminger and Labrador Seas, ASTE_R1 can reliably reproduce the observed hydrographic variability (Figures 14a and 14b). However, the solution exhibits a widespread systematic warm bias between 500 and 2,000 m underlying a cold bias in the upper ∼500 m over the subpolar region (shown for the Labrador and Irminger Seas in Figure 14c). Overall, salinity is biased fresh in the Labrador Sea, whereas a salty bias characterizes the Irminger Sea ( Figure 14c) and the wider eastern subpolar North Atlantic region (not shown).
In the subtropical North Atlantic, a large fraction of the salinity misfit in it0 is due to an excess freshwater flux from the atmosphere. Comparison to the independent Global Precipitation climatology Project version 2.3 product (GPCPv2.3, Adler et al., 2018) reveals an excess precipitation bias in JRA55, that is most pronounced in the North Atlantic and subpolar gyre region of the ASTE domain (Figure 13d), and that resulted in a large fresh bias in the upper ∼500 m of the unconstrained it0 solution (Figure 13c). The adjoint-based optimization provided a systematic approach for removing this excess precipitation bias, such that after ∼12 iterations the misfits to Argo salinity in the upper ocean reduced to within the observed uncertainty NGUYEN ET AL.  Figure 13b). Consequently, this improvement also yields better agreement in the mean with an independent GPCPv2.3 data set ( Figure 13d). It is important to stress again that these adjustments are made whilst retaining the ocean model dynamical and kinematical consistency.
The published literature offers notable differences in tracer reference values employed in the computation of reported heat, freshwater and volumetric watermass transports. These range from regional basin means (Beszczynska-Möller et al., 2012;de Steur et al., 2018;Smedsrud et al., 2010;Tesdal & Haine, 2020) to gateway and surface means (Tsubouchi et al., 2018) to freezing temperature in the Arctic (Beszczynska-Möller et al., 2012;Woodgate, 2018). In some cases, transports were computed along a particular range of isopycnals (Tsubouchi et al., 2018). Heat transport computed in ASTE_R1 assumes a reference temperature θ r = 0°C (most accurate numerically, see Appendix B). For the Bering Strait, we also compute heat transport referenced to the freezing temperature of seawater θ r = −1.9°C to facilitate comparison with published estimates. For the computation of freshwater transports, we assume a reference salinity S r = 34.8 ppt and integrate from the surface down to the reference isohaline. We refer the reader to Appendix B for details on potential errors incurred when computing transports using non-zero reference values. Due to the difference in reference values employed here and in some of the studies listed above, we seek consistency in terms of comparable transport magnitudes as opposed to exact agreement.
To provide a useful comparison of ASTE_R1 mean transports with those reported in the literature, it is important to note whether published estimates are based on historic data or more recent acquisitions, given how fast the high latitudes are observed to be changing. In the first four years of the ASTE_R1 period  Note. All uncertainties provided are given in terms of standard deviations based on monthly estimates after the seasonal climatology has been removed. FW transport is computed using Equation B3.2 of Appendix B. a The vertical convergence of FW, from air-ice-sea fluxes, is the same as that for volume and is exact. Lateral convergence and tendency of FW, however, are approximate. As a result, the budget for FW is not fully closed (see Appendix B). b Franz Josef Land, c Severnaya Zemlya.  -Möller et al., , 2012Curry et al., 2014;de Steur et al., 2014;Hansen et al., 2015;Skagseth et al., 2008;Woodgate, 2018), and modeling studies (Heuzé & Årthun, 2019;Ilicak et al., 2016;Tsubouchi et al., 2018). In addition to net transports, we also provide estimates of transports of important watermasses at Fram Strait (as defined in Beszczynska-Möller et al., 2011) and through the GSR (as defined in Hansen & Østerhus, 2000;Østerhus et al., 2019). For the in/outflow transport estimates given in Figure 15 Figure 15), southward surface flow of freshwater across DS and dense overflow across the entire ridge (green and blue color in Figure 15).
Watermass definitions for surface outflow, dense outflow, modified water, and inflow AW in ASTE_R1 can differ from Østerhus et al. (2019) and Hansen and Østerhus (2000) for the reasons outlined in Section 2.3. Our choice for σ θ is justified in Appendix C. For the overflow through DS, the range of 27.4 ≤ σ θ ≤ 27.8 used in ASTE_R1 is associated with southward transports of −1.6 ± 0.9 to −0.5 ± 0.3 Sv (shown in blue in Figure 15), corresponding to 16%-50% of the observed estimate using σ θ = 27.

Reference (T r , Period)
Gateway the observed estimate at the DS by ∼30%. In total, the net volume transport across DS in ASTE_R1 is about 47% of that reported by Østerhus et al. (2019).

For the Arctic Ocean and GIN Seas heat and freshwater budgets, transports through Fram Strait (FS) and
the Barents Sea Opening (BSO) are also important. Across FS, the inflow of warm AW along the West Spitsbergen Current (red color in Figure 15) is 6 ± 1 Sv in ASTE_R1, with 3.4 ± 1.1 Sv carrying the core AW water warmer than 2°C. This is consistent with corresponding estimates of 6.6 ± 0.4 and 3.0 ± 0.

Heat Transports
All ASTE_R1 net heat transports are northward into the Arctic Basin. Time-mean transports across key gateways are consistent with observation-based estimates, a result that is aided-although by no means guaranteed-by constraining the state estimate using mooring T/S data (Table 4, Figure 16).
At Bering Strait ASTE_R1 heat transport is 14 ± 8 TW (referenced to T r = −1.9°C), consistent with the 11.6 TW to 14.3 TW range determined by Woodgate (2018). At Davis Strait, the ASTE_R1 estimate of 20 ± 4 TW is consistent with 20 ± 9 TW obtained by Curry et al. (2011). Across the GSR, heat transport is in good agreement with previous published estimates across the two eastern channels (IF and FS, Figure 16), but is underestimated across the Denmark Strait. Here, the total poleward diffusive heat flux dominates and opposes the equatorward advective term in ASTE_R1. This diffusive dominance has also been suggested using heat budget analyses in ECCOv4 (Buckley et al., 2015). A more detailed discussion of the full time-series and contributions of advective and diffusive fluxes to the total transport is given in Appendix B.
Further north, at Fram Strait, ASTE_R1 poleward heat transport is 39 ± 12 TW (referenced to T r = 0°C), consistent with 36 ± 6 TW from Schauer and Beszczynska-Möller (2009)  The Fram Strait heat transport is increased by approximately one third (15 TW) on accounting for sea ice advection. Heat transport into the Barents Sea across BSO of 62 ± 15 TW is consistent with observation-based estimates of between 48 TW and 73 TW (Rossby et al., 2018;Skagseth et al., 2008;Smedsrud et al., 2010). Most of this heat is lost via air-sea exchange in the Barents and Kara Seas (Lind et al., 2018), yielding negligible heat transports from this shallow region into the Arctic Basin ( Figure 16). Air-sea exchange also accounts for significant loss of heat in the Canadian Arctic Archipelago, such that only ∼35% of the amount transported across the Davis Strait reaches the Arctic Basin.
There is a large spread amongst existing observation-and model-based studies with significant disagreements even after accounting for uncertainty (Figure 16), due in part to the lack of common data period and reference temperature used in the calculations. Overall, nevertheless, the poleward heat transports in ASTE_R1 are in good agreement with previous estimates (Figure 16).

Freshwater Transports
The practice of reporting ocean freshwater (FW) transport/content in place of absolute salt transport/content is ubiquitous in the literature, but plagued by the need to specify a reference salinity, S r , and to choose the vertical extent over which the integral is computed (i.e., full depth vs. to the depth of the reference salinity, S r z ). No unique choice emerges from consideration of seawater physics. Instead, S r is selected inconsistently between studies. As cautioned by Schauer and Losch (2019), this not only complicates comparisons but can gives very different impressions of the changing ocean state, due to strong sensitivity to the choice of S r . Acknowledging this issue, we nevertheless elect to report FW transport below (as did Tesdal and Haine (2020) in their recent study focusing on the subpolar North Atlantic and Nordic Seas), in order to conduct our assessment of ASTE_R1 hydrography in the context of existing estimates. To the best of our knowledge no published observational estimates report Arctic salt transports that would provide a basis for comparison (the modeling study by Treguier et al., 2014 is one known exception). We proceed with caution and flag comparisons for which calculations differ. ASTE_R1 FW fluxes are reported using a reference salinity of S r = 34.8 ppt and integrated down to S r z . Our calculation uses monthly averages of both the Eulerian velocity and salinity. In Appendix B, we provide a detailed discussion of the potential errors in FW calculations with these choices, along with errors incurred in omitting bolus and diffusive terms. Salt transport or salt content changes in ASTE_R1 will be revisited in future work.
Similar to our assessment of volume and heat transports, we start by examining FW transports across the gates into the Arctic Mediterranean. At the Bering Strait, ASTE_R1 combined liquid and solid FW import of 1,711 ± 608 km 3 /yr is lower than the 2,670 ± 144 km 3 /yr estimated by Woodgate et al. (2015) and Woodgate (2018). There are several candidates to explain this ∼960 km 3 /yr FW transport deficit at this gate. Most importantly, the river runoff climatology used in ASTE_R1 has likely not taken into account potential increased discharge from the Yukon River into the Bering Sea just upstream of the strait ( Spreen et al. (2009), respectively. A main reason for the lower (by ∼1800 km 3 /yr) liquid plus solid export in ASTE_R1 across Fram Strait is the lower (by ∼960 km 3 /yr) net FW import through the Bering Strait relative to observations. An additional inconsistency is the use of a runoff climatology in ASTE_R1, which fails to account for Greenland solid/liquid discharge and its observed recent increase into the Arctic sector by ∼105 km 3 /yr (Bamber et al., 2012), as well as increased river outputs (as reported in Bamber et al., 2012;Proshutinsky et al., 2020). With respect to the latter, ASTE_R1 has a deficit of ∼220 km 3 /yr. The climatology also does not account for Greenland FW (combined solid and liquid) discharge into the GIN Seas and Baffin Bay of nearly 150 and 250 km 3 /yr, respectively (Bamber et al., 2012).
In the Canadian Arctic Archipelago, ASTE_R1 has a FW flux deficit of ∼226 km 3 /yr from land ice (Carmack et al., 2016). These omissions likely contribute to the underestimation of southward FW transports across both the Denmark Strait (by ∼940 km 3 /yr) and Davis Strait (by ∼150 km 3 /yr).
The net lateral convergence of the combined liquid and solid freshwater flux in ASTE_R1 of −2,510 ± 1272 km 3 /yr is nearly balanced by the net vertical convergence of 2,353 ± 324 km 3 /yr, yielding a net tendency of −298 ± 1,131 km 3 /yr. For the liquid flux alone, the tendencies are −1,484 ± 1123 for lateral convergence, 1,238 ± 2,478 vertical convergence, and −284 ± 2,169 km 3 /yr total tendency.

Heat and Freshwater Storage
Complementing the transport estimates, we conclude our initial assessment of ASTE_R1 with an overview of derived basin-scale 2002-2017 time-mean and time-variable heat and freshwater budgets, focusing on comparisons between ASTE_R1 and existing estimates. A full assessment of the mechanisms underlying Arctic Mediterranean and subpolar North Atlantic heat and freshwater content change over the ASTE_R1 period will be addressed in a separate study. As noted earlier, decisions made in our tracer transport/budget calculations facilitate these comparisons but are non-unique. For freshwater transports/budgets this introduces ambiguity that is best resolved prior to detailed dynamical investigation.

Heat Content
Considering the Arctic region of ASTE_R1 in its entirety, the net heat input from convergence of horizontal (ocean plus ice) heat transports (53.3 ± 12.0 TW) exceeds the net heat loss due to local air-ice-sea fluxes (−27.4 ± 28.7 TW) by a factor ∼2, yielding a net heating rate of 25.9 ± 30.0 TW when averaged over the period 2006-2017 (Table 4). Relative to horizontal convergence, vertical exchange at the air-ice-sea interface is significantly less well constrained due to large uncertainties in atmospheric reanalyses at high northern latitudes (Beesley et al., 2000;Chaudhuri et al., 2014;. Partitioned by basins, similar enthalpy gains are estimated for both the western (14.0 ± 24.3 TW) and eastern (12.5 ± 14.0 TW) Arctic region. In the water column, heat gain is concentrated mainly in the AW layer (240-1,000 m) in both the western (9.3 ± 2.7 TW) and eastern (6.4 ± 4.6 TW) Arctic. In the upper 60 m of the water column, the tendency is negligible but with large variability (0.2 ± 25.3 TW) due to mixed layer processes and exchange with the atmosphere.
NGUYEN ET AL.

10.1029/2020MS002398
26 of 49 Warming since the early 2000s has been reported in the Arctic, documented alongside enhanced "Atlantification" in the Eastern Arctic (Polyakov et al., 2017(Polyakov et al., , 2020) and a fivefold increase in solar absorption by near surface waters in the Western Arctic (Jackson et al., 2010;Timmermans et al., 2018). This warming proceeds at a sustained rate, with recent studies suggesting a doubling of OHC in the Beaufort Gyre halocline between 2003 and 2013 . Figure 18a shows a comparison of the time-series of the halocline heat content in ASTE_R1 and estimates based on ITP and mooring data in the Beaufort Sea. As discussed in Section 2.3, adopting exact watermass classifications from observational studies may be inappropriate for model analysis, due to representation error of subgrid scales. Thus in addition to the salinity limits used to identify the upper halocline layer of 31.0 ≤ S ≤ 33.0 in Timmermans et al. (2018), we also compute the heat content sensitivity to the salinity bounds. By changing the near-surface lower salinity bound within the range 31-31.5 ppt, the mean halocline heat content changes by 2%-3% per 0.1 ppt increment, but the variability and trend remain unchanged. This confirms that despite the systematic warm bias, the positive trend in halocline heat content is well-captured in the ASTE_R1 solution.
In the Canada Basin, ASTE_R1 exhibits a warming rate of 9.3 ± 2.7 TW for the period 2006-2017 (Figure 16). There are insufficient observations to validate this directly, but we can corroborate our estimate with a back of the envelope calculation as follows. Recent ITP acquisitions report core AW temperatures in this region ∼0.5°C warmer than the PHC climatology (Figure 19a), where the latter is representative of the second half of the 20th century. Since ITPs only measure to ∼800 m depth, we conservatively assume a depth-average warming of 0.20-0.25°C over the 170-1,000 m range in the Western Arctic basin interior (area 4400 × 10 3 km 2 ), yielding a warming rate of ∼5.7-7.2 TW, about two thirds of the rate estimated in ASTE_R1. Compared to ITP data, ASTE_R1 shows a positive bias of ∼0.15°C in the core AW temperature NGUYEN ET AL. in the Western Arctic basin, which accounts for the higher tendency here. However, we note that this ASTE_R1 bias is within the combined data and representative error σ T , which is not the case for the PHC bias ( Figure 19b).
The Eastern Arctic suffers from an extreme paucity of data, such that even back of the envelope estimates of basin-wide heat content (and its tendency) are not possible. Instead, we turn to recent observations for evidence of warming in this basin. Pulsed injection of AW at Fram Strait has been documented by Polyakov et al. (2011). A notable warm anomaly pulse of ∼1°C entered the Arctic in 2004. It has subsequently been observed crossing the NABOS section at 126°E, and has been recorded further downstream at numerous sections along the eastern basin's rim (Polyakov et al., 2011). In addition, the seasonal amplitude within the halocline has been observed to increase by 0.75°C between 2004 and 2015 (Baumann et al., 2018;Polyakov et al., 2017). These observations provide evidence for a warming Eastern Arctic and a weakened halocline. The latter is accompanied by shoaling of the AW layer toward the bottom of the mixed layer and increasing heat ventilation (Polyakov et al., 2020). This is a mechanism by which heat along the AW pathway is removed instead of being sequestered at depth. To determine the relative importance of these two mechanisms (ventilation vs. sequestration) in contributing to the positive heat content tendency at different depths and throughout the Arctic in ASTE_R1, a more detailed analysis of AW circulation and ventilation will be needed in future work.
In the Barents Sea, Lind et al. (2018) documented pronounced increases in decadal mean OHC in the upper 100 m of the water column, which they attributed to an increase in AW inflow through the Barents Sea Opening. In Figure 18c

Freshwater Content
Based on observations, predominantly from satellite altimetry and ITPs in the Beaufort Sea, the liquid freshwater content (FWC) in the Arctic has been estimated to be increasing . Proshutinsky et al. (2020) summarized recent works attributing this FWC increase to several factors, including shifts in atmospheric circulation, increased FW fluxes through the Bering Strait, and increased runoff from the MacKenzie river. A comparison between ASTE_R1 and Proshutinsky et al. (2019) estimates for FWC in the Beaufort Gyre shows that ASTE_R1 captures the observed increase in FWC between the 2004 and 2008. With the exception of a small decrease in 2015-also seen in the observations-the Beaufort Gyre FWC in ASTE_R1 remains relatively constant for the period 2008-2017 (Figure 18b). Proshutinsky et al. (2019) report an increase from 2015 to 2017 which is likely missing from ASTE_R1 due to the omission of both increased river runoff and land ice discharge in our forcing climatology and absence of the observed increase in FW import through the Bering Strait as previously discussed (Figure 17).
The connection between the Arctic FWC increase and circulation changes in the GIN Seas and North Atlantic has been the subject of several investigations (Carmack et al., 2016;Dukhovskoy et al., 2016;Tesdal & Haine, 2020). A recent review by Haine et al. (2015) summarized Arctic exchanges with the Canadian Arctic Archipelago (north of Davis Strait) and the Barents Sea (east of the Barents Sea Opening). These estimates heavily rely upon atmospheric reanalyses (for the provision of surface fluxes) and unconstrained model output (for tracer content change). Further south, Dukhovskoy et al. (2019)    GIN Seas, highlighting large uncertainty due both to lack of constraint and acute dependency of transports on model resolution (see also Weijer et al., 2012).

Discussion
A preliminary discussion of how the optimization acts to bring the model into consistency with the available observations focuses on the question of which control variables played a dominant role in achieving a reduction in misfit (Section 5.1). A second point of discussion highlights known issues with this first release of ASTE and suggestions on how to improve future releases (Section 5.2).

Identifying Key Control Adjustments
Production of the ASTE_R1 solution was achieved by gradient-based optimization, which iteratively adjusts a set of control variables (Section 2.2). Our control space (Ω) comprises 3D fields of initial (i.e., January 1, 2002) temperature and salinity (θ 0 , S 0 ), time-mean spatially varying ocean mixing coefficients (  , gm   , and z  ), and time-varying 2D fields of near-surface atmospheric state variables (T air , q air , u w , v w , R sw , R lw , and P).
We now seek to identify which of these control variable adjustments had the largest impact on reducing the model-data misfit in ASTE_R1 relative to the unoptimized it0.
To proceed, we performed two forward ensemble experiments, each experiment consisting of eight members. In the first experiment, individual optimized control variables from ASTE_R1 were substituted into it0, which was then re-run. Each ensemble member is characterized by containing one of the ASTE_R1 optimized control variables or variable pairs listed in Table 5 (left two columns). Note that there are 8 variables that come in pairs, thus making 4 total pairs: the optimized initial conditions ( However, a note of caution is needed. The control variables are not fully independent (e.g., some of the atmospheric state variable controls are related via bulk formulae or shared physics), and as a result, it is not possible to determine their full impact in isolation. For this reason, we performed a second set of experiments, reversing the sense of the substitutions, such that the non-optimized controls from it0 were substituted into ASTE_R1, which was then re-run (i.e., the optimized control variables were reset to their first guesses). In this experiment, each ensemble member is characterized by containing one of the it0 non-optimized control variables or variable pairs listed in Table 5 (right two columns). This second experiment lends confidence to our assessment as follows: an optimized control is highly likely to be an important ingredient of the ASTE_R1 solution if its incorporation notably improves the it0 re-run (first ensemble) while its omission notably degrades the ASTE_R1 re-run (second ensemble).
In Figure 20 we examine the impact of the control substitutions on the costs in both ensemble experiments. We show only normalized costs with respect to the following aggregated data sets: Argo, ITP, Beaufort Gyre moorings, and satellite-based observations of SST, SSH and sea ice concentration. This choice enables a more focused discussion whilst also informing the large-scale quality of the solution near the ocean surface in the Atlantic Ocean (where the majority of SSH and SST data were acquired) and throughout the upper ocean in the North Atlantic, GIN Seas, and Labrador Sea interior (from Argo T and S data) and in the western Arctic (from ITP and Beaufort Gyre moorings). Lastly, costs for sea ice concentration indicate performance of modeled air-sea fluxes and mixed layer properties in marginal ice zones (see Figure 7 and related discussion).
NGUYEN ET AL.
Note. In Experiment 1, the control variables listed in column #2 were added to unoptimized it0; in Experiment 2, the control variables listed in column #4 were withheld from optimized ASTE_R1.

Table 5 Ensemble Members for Each of the Two Ensemble Experiments
Our analysis, based on Figure 20, reveals the importance of precipitation (P) adjustments in obtaining realistic subsurface salinity distributions in the North Atlantic through the removal of the systematic excess rain bias discussed in Section 3.3 ( Figure 13d). Specifically, inclusion of the optimized precipitation P R1 in the it0 re-run (column under "P" on the left half of Figure 20) reduced the cost with respect to Argo salinity by 56% (orange square at row "Argo-S″ and column "P" corresponding to large negative, that is, cost reduction, values between −50 and −60 as indicated in the color scale). A 31% increase in this cost was seen on omission of the optimized precipitation in the ASTE_R1 re-run (blue square at row "Argo-S″ and column "P" on the right half of Figure 20, corresponding to positive, i.e., cost increase, values between +30 and +40 in color scale). Large improvements in SSH can also be attributed in part to corrected precipitation, as well as surface winds. Amongst other atmospheric forcing variables, surface air temperature, downward radiative forcing, and winds all have important impact on the sea ice cover (collective average of 17% improvement to it0 and 5% degradation to ASTE_R1).
Adjustments to the initial conditions, vertical diffusivity and eddy mixing are also found to be important for improving subsurface hydrography throughout the ASTE_R1 domain. Optimized eddy mixing-related controls (  1 1 , R R gm   ) alone result in a 14% improvement (and 14% degradation) of the Arctic hydrography misfit when included (omitted) from the it0 (ASTE_R1) re-runs, respectively. These adjustments to the eddy mixing parameters also improve SST and SSH, and-in addition to the adjustments made to the vertical diffusivity and atmospheric conditions-are also seen to be critical for improved representation of sea ice cover.

Known Issues and Future Directions
During production of ASTE_R1 we have striven to utilize all constraints known to us and that the state estimation machinery could handle. This comprises O(10 9 ) observations from diverse data sources (Table 2).
NGUYEN ET AL.

10.1029/2020MS002398
30 of 49 Figure 20. Percentage change (color) in cost with respect to the constraint listed on the ordinate, attributable to the control substitutions given on the abscissa. The left group are ensemble members of perturbation Experiment 1, for which optimized controls Ω R1 are substituted into the it0 re-runs. The right group are ensemble members of Experiment 2, for which non-optimized controls Ω i0 are substituted into the ASTE_R1 re-runs. On the left half of the plot, negative values indicate an improved solution i.e., a cost reduction with respect to it0, implying that the optimized controls are important for reducing the misfit. On the right half, positive values indicate deterioration of the solution i.e., a cost increase with respect to ASTE_R1, implying that these control adjustments are critical for obtaining the optimized ASTE_R1 state. Colored rectangular outlines highlight patterns of most impactful control variables (e.g., precipitation is important for the reduction of costs to Argo salinity and SSH). Despite this effort, some systematic biases remain in the ASTE_R1 solution. As the optimization is ongoing and ASTE is still converging, we anticipate the costs listed in Table 3 will continue to reduce and some of the remaining biases will be removed. In certain cases, due to model structural errors or non-resolved physics, full convergence might not be attainable (Wunsch & Heimbach, 2007). Here we discuss notable issues remaining in the ASTE_R1 solution and possible future directions for developing the next ASTE release with improved model physics.
Eastern Arctic hydrography: One of the largest remaining systematic biases is found in the Eurasian Basin, where subsurface constraints comprise sparse ITP sampling of the basin interior alone. Although the inflow is constrained by moorings at Fram Strait, downstream observation of the circulation, eddy-induced stirring and vertical mixing in the Eurasian Basin along the shelf-basin slope and interior are limited (Figure 21a). The serious implications of this paucity of data are highlighted by considering that the AW inflow takes ∼6-10 years to transit this region, during which there are no local observational constraints. As a result, the inverse problem is highly under-determined. In practice, under-determination allows non-unique pathways to misfit minimization. For ASTE_R1, we find that the AW layer in the Eastern Arctic spreads to occupy a greater depth range toward the end of the estimation period. This problem was also present in it0 and has been partly ameliorated during the optimization, as reflected in the removal of the largest positive temperature misfits with respect to ITP data (Figures 21b and 21c). This thickening of the AW layer is also a common problem in many state-of-the-art Arctic Ocean models (Docquier et al., 2019;Holloway et al., 2007;Ilicak et al., 2016;Uotila et al., 2019).
As it is unlikely that widespread observation of 3-D velocity and mixing will be made in the foreseeable future, we anticipate that AW watermass representation in the Eurasian basin will remain an issue for both the next generation of state-of-the-art Arctic Ocean models and the next ASTE release. Although we do not expect large gains from planned changes to the ocean observing system in the near future, we do anticipate improvements in sea ice state, mixed layer representation, and shelf-basin exchanges in the next ASTE release, due to recent improvements to the stability of the adjoint of the sea ice thermodynamics (Bigdeli et al., 2020). This will enable a more complete use of sea ice observations as active contributions to the cost function reduction (Equation 1). The sensitivity of the associated model-data misfits to the control space can then be used to better adjust atmospheric forcings. This will allow us to fully leverage the constraint from satellite-based observations of the sea-ice state, which could only be partly exploited in the pseudo sea ice adjoint employed for construction of ASTE_R1. Inclusion of the sea ice thermodynamics adjoint could NGUYEN ET AL. potentially improve AW upward ventilation (Ivanov et al., 2012;Polyakov et al., 2020) and preserve a more stable AW layer thickness, both in the Eurasian Basin and further downstream in the Western Arctic.
Arctic Circumpolar Current: In the Laptev Sea it is thought that the circumpolar circulation of AW splits at ∼145°E, with a fraction returning to Fram Strait along the Lomonosov Ridge (Rudels, 2015) and the remainder continuing along the basin's rim into the Western Arctic, although the exact partitioning is not well constrained. In the Western Arctic it is typically assumed that the AW continues to circulate cyclonically along the basin boundary, although both ASTE_R1 (Grabon, 2020) and a modeling effort informed by observed radionuclide distributions (Karcher et al., 2012) suggest a weak anticyclonic circulation during the last decade. Recent work analyzing all available current meters (updated from Baumann et al., 2018) has yielded velocity probability distributions for the Arctic region. This will be investigated as a novel approach to constrain ocean velocities within the AW circulation in the next ASTE release. In addition, a more detailed examination of the momentum and vorticity budgets along the circumpolar current will offer insights into the role of viscous dissipation and eddies in maintaining the cyclonic sense of circulation (Spall, 2020;Yang, 2005).
Arctic river runoff and Greenland discharge: FW transports and content in the late 2010s are low in ASTE_R1 relative to independent observations (Section 4.3 and 4.4.2). Near the surface in the Arctic and along the Greenland coast, recent increases in river (Shiklomanov et al., 2020) and tundra runoff (Bamber et al., 2012), surface solid and subsurface glacial discharge (Bamber et al., 2012(Bamber et al., , 2018 have been observed. This increase was not included in the ASTE_R1 forcing. Meaningful application of these FW fluxes as model forcings, especially in the Arctic marginal seas, requires careful consideration of the following factors. Sub-glacial discharge is observed to enter the outlet glacier fjord at depths near the grounding line instead of at the surface of the fjord's exit to the continental shelf (Sciascia et al., 2013;Straneo & Cenedese, 2015). Mixing and entrainment of this FW with the surroundings creates modified water whose property is prohibitively difficult to continuously track downstream from the source using observed T/S (Beaird et al., 2018). Consequently, the pathways of FW redistribution are highly uncertain. Numerical simulations with Greenland discharge distributed at the surface yield pathways from the source into the interior of SPNA and GIN seas that vary substantially with model resolution and representation of mean currents (Dukhovskoy et al., 2016;Weijer et al., 2012). The depth to which this FW is mixed down also varies highly with resolution (Dukhovskoy et al., 2016), causing near surface over-freshening in certain cases and a 30%-50% decrease in the North Atlantic Meridional Overturning Circulation (AMOC) at time-scales varying between 3 and 50 years (Weijer et al., 2012). Similarly, preliminary sensitivity experiments in ASTE_R1 with observed Greenland discharge applied at the surface show over-freshening of the upper ocean in the Greenland Sea and a decrease in the AMOC at 55°N by 40% within 5 years, inconsistent with observations (not shown). Prior to the next ASTE release, a dedicated study will be required to implement updated estimates of Greenland discharge as a subsurface freshwater forcing, consistent with observations (Straneo & Cenedese, 2015). This will entail incorporation of a melt water plume parameterization into the ASTE framework. Lastly, instead of being absorbed into net surface freshwater flux E − P − R, a new control variable for runoff could be introduced to isolate and fully interrogate sensitivity to subsurface forcing from subglacial discharge.

Subpolar North Atlantic hydrography:
A warm bias in ASTE_R1 at 500-2,000 m depth persists both in the Irminger Sea (Figure 14c) and throughout the eastern SPNA (not shown), and is associated with a weaker poleward transport of Atlantic warm water across the GSR (Figure 15). Poor representation of AW inflow across the GSR is a common problem in coarse to medium resolution ocean models (Heuzé & Årthun, 2019). Specifically, these models produce lower volume and heat transports across the GSR compared to observations (Heuzé & Årthun, 2019). Since the resolution of ASTE_R1 is ∼18 km in the subpolar gyre, we anticipate incomplete representation of both eddy/diffusive mechanisms-estimated to be important across the shallow Denmark Strait and Iceland-Faroe Ridge (Buckley et al., 2015)-and watermass transformations in the ASTE_R1 solution. Thus, although ASTE_R1 can capture the mean transports of volume and heat between Iceland and Scotland (Figures 15-17), the remaining warm bias across Denmark Strait and south of the GSR likely impacts our estimate of heat content in both the eastern SPNA and Nordic Seas and alters the optimized air-sea heat flux in both regions. Recent data from the Overturning in the Subpolar North Atlantic Program (OSNAP) observing system (Lozier et al., 2017(Lozier et al., , 2019  will provide important information in the subpolar region, and an especially valuable constraint on the boundary currents and overflow waters, not captured by Argo.

Summary and Outlook
We have presented the first release of the Arctic Subpolar gyre sTate Estimate, ASTE_R1, a data-constrained and dynamically consistent ocean-sea ice synthesis spanning the period 2002-2017. ASTE_R1 is produced using the ECCO adjoint-based state estimation framework, in which an ocean general circulation model, the MITgcm serves as a dynamical interpolator, spreading the influence of O(10 9 ) incorporated observations through space and time by way of linearized adjustment processes encapsulated in an adjoint model. Importantly, the model-data misfit is reduced via iterative adjustments to the initial hydrographic conditions, atmospheric forcing and model mixing parameters alone, ensuring adherence to the governing equations throughout the entire estimation period. This distinguishes our approach from ocean reanalysis, in which violation of conservation laws complicates application for climate research (Stammer et al., 2016). The ability to assess closed tracer and momentum budgets in ASTE_R1 is a key strength of the product. As all sources and sinks are accounted for, full heat, salt and momentum (or vorticity) budgets can be analyzed to identify dominant sources contributing to the observed changes. These closed budget analyses can also be performed in T, S, and σ space following Abernathey et al. (2016), enabling diagnosis of watermass evolution and destruction in the ASTE_R1 solution. In addition, the adjoint modeling infrastructure allows for linear sensitivity studies using ASTE_R1 for investigation of causal mechanisms underlying variability in key quantities of climate interest (e.g., Bigdeli et al., 2020;Nguyen, Woodgate, & Heimbach et al., 2020;Pillar et al., 2016).
During production of ASTE_R1 we have strived to utilize all observational constraints known to us and that the state estimation machinery can handle. ASTE_R1 thus arguably represents the biggest effort undertaken to date with the aim of producing a specialized Arctic ocean-ice estimate, freely available to the research community. This complements existing global ECCO solutions (Forget, Campin et al., 2015;, the Southern Ocean State Estimation (SOSE, Mazloff et al., 2010) and other global and regional ECCO derivatives (e.g., Gopalakrishnan et al., 2013;Köhl, 2020;Köhl & Stammer, 2008;Zaba et al., 2018).
For this initial assessment of ASTE_R1, we have focused on comparison to available observational constraints. Many of these were actively employed in the optimization procedure, but some (e.g., all volume and tracer transport estimates) were withheld, allowing independent verification. The optimized solution serves as a significant improvement from the unconstrained state, achieving consistency with the majority of incorporated observations, including both the set used in the optimization and that retained for post-validation (Table 3).
The most substantial misfit reduction in ASTE_R1 are sea ice cover in the marginal ice zone, western Arctic hydrography, and subtropical North Atlantic sea level anomaly and subsurface salinity (Table 3, Figure 6). In the Arctic Mediterranean, using only a proxy sea ice adjoint, ASTE_R1 achieves a 83% reduction in misfit to satellite-derived sea ice concentration constraints, mainly via improved representation of the sea-ice edge (Figure 7). The solution faithfully reproduces both the observed seasonal cycle and low frequency trend of sea-ice extent.
At Fram Strait, the mooring array is crucial to constraining the important AW inflow and local hydrographic properties. At this important Arctic gateway, ASTE_R1 exhibits a 58% misfit reduction through the water column across the strait relative to the unconstrained simulation. In the Arctic interior, ITPs provide unique information on the subsurface hydrography. Because 71% of the ITP profiles are located within the upper 5-800 m in the Canada Basin interior, the most significant misfit reduction was seen here (85% in T and 62% in S, Figure 10). In the remaining Arctic basin, low data coverage, combined with large uncertainty in the mean circulation and mixing parameters, resulted in less notable improvement (reductions of 89% in T and 31% in S), but biases persist, especially at depth below the AW core ( Figure 21).
Accompanying improved fit to hydrographic data used to constrain the solution, we find improvements in basin-scale heat and freshwater content. Interannual variability and low frequency trends in both heat and FW content are well represented in the Arctic Mediterranean and SPNA of ASTE_R1. In the Beaufort Sea, NGUYEN ET AL. We have been careful to clearly outline the notable biases remaining in the ASTE_R1 solution. These include a warm bias below the AW core in the eastern Arctic and in the east SPNA. The cause is a combination of lack of constraint here for the hydrography, mean circulation, and the adjustable initial condition and parametric controls. Additional biases exsists in FW transports and contents in the Arctic Mediterranean due to the omission of increased runoff from Arctic rivers and Greenland freshwater discharge.
An advantage of our approach is that the use of a dynamical interpolator can improve spectral representation of the estimated state compared to gridded products produced using statistical interpolations (e.g., Verdy et al., 2017). This has not been addressed here, but it is a useful avenue for future ASTE assessments and ongoing development.
Looking toward the next ASTE_R1 release, we expect the greatest progress will be made by incorporating new model physics. In particular, improving the stability of the sea ice thermodynamic adjoint (Bigdeli et al., 2020) will enable its use in ASTE, providing stronger constraint of air-ice-sea exchanges and ocean ventilation. Future efforts will target hydrographic improvements along the Arctic shelf-basin slope in the eastern Arctic to reduce the ASTE_R1 AW layer warm bias. Additionally, updated estimates of runoff and calving fluxes and inclusion of a parameterization of sub-glacial discharge will enable improved estimate of freshwater redistribution and interbasin exchange. New constraints, including datasets from the OSNAP mooring array (Lozier et al., 2017(Lozier et al., , 2019, sea surface salinity (Fournier et al., 2019;Vinogradova & Ponte, 2012), and sea ice thickness (Ricker et al., 2017;Tian-Kunze et al., 2014) will also be fully utilized in the production of a further improved next ASTE release.

A1. Configuration Set Up
The model configuration and all necessary inputs, including the optimized control adjustments, required for ASTE_R1 re-runs are available to the public, as discussed in the next section. The code base employed for ASTE_R1 production was MITgcm checkpoint c65q. ASTE_R1 was built using the full state estimation infrastructure, including specialized packages for misfit and adjustment evaluation, developed for ECCOv4r1 (Forget, Campin et al., 2015). In addition, two code developments specific to ASTE include the implementation of a vertical diffusivity power control ( 10 log z  ) and the capability to switch between daily and monthly SSH costs.
To ensure numerical stability during ASTE_R1 production, the following model choices were important: (a) a staggered time-step for momentum advection and Coriolis terms; (b) third-order advection for tracers (scheme code 30 in Table 2.2 in Adcroft et al., 2018), (c) linear free surface approximation, and (d) application of freshwater forcing via a virtual salt flux (i.e., no accompanying change in mass). These choices permitted a time-step of 1200 s. After 62 iterations, better model choices were used for the final forward run that produces more accurate physics in the distributed version of ASTE_R1. These include (i) seventh order advection for tracers (scheme code 7 in Table 2.2 in Adcroft et al., 2018); (ii) nonlinear free surface with scaled z* coordinates , and (iii) application of freshwater forcing via a real freshwater flux (i.e., with accompanying change in mass, Campin et al., 2004). These choices required a shorter time-step of 600 s. The ASTE_R1 solution described and assessed in this paper is from the re-run of iteration 62 with the model choices (i)-(iii) described above.
In the distributed code, at compile and run-time, the user has the choice to use the more stable set up with a time-step of 1200 s or employ the more accurate numerics and physics with a time-step of 600 s as described above. We found that these small changes in the model configuration for the final forward run did not have a significant impact on the solution. This result is consistent with published studies suggesting small differences in ocean dynamics between LFS versus NLFS in combination with virtual salt or real freshwater flux Roullet & Madec, 2000;Yin et al., 2010). The advantage of their application here is in enabling more accurate physical interpretation of mass and freshwater budgets. However, since these options also require a shorter timestep (for the nonlinear free surface) and a larger stencil (for the higher NGUYEN ET AL.

Journal of Advances in Modeling Earth Systems
order advection), their use demands significantly more computational resources (twice the wallclock time). For this reason, it was not feasible to employ these options until the final stages of ASTE_R1 development.

A2. Distribution of the ASTE_R1 Solution
The full ASTE_R1 solution is publicly available through the NSF Arctic Data Center as follows: a) Time varying fields as monthly averages and monthly snapshots (Nguyen et al., 2021a); b) Depth-integrated time varying fields as monthly averages and monthly snapshots (Nguyen et al., 2021b); c) Selected time varying state variables as daily averages (Nguyen et al., 2021c); d) 12-months climatological averages (Nguyen et al., 2021d); e) In situ profiles and model-equivalent (Nguyen et al., 2021e); f) ASTE_R1 Grid files, Documentations (user guide, domain layout) and MATLAB toolbox to help analyze the output fields (Nguyen et al., 2021f); g) Compile time and run time inputs necessary to reproduce ASTE_R1 with the MITgcm (Nguyen et al., 2021g).
All model output fields are available here as NetCDF files. In addition to being archived at the Arctic Data Center, ASTE_R1 NetCDF data are also mirrored at the UT-Austin ECCO portal at: https://web.corral.tacc. utexas.edu/OceanProjects/ASTE/, which is provided by the Texas Advanced Computing Center (TACC). Alternative to NetCDF format, the monthly mean fields are additionally hosted in a compressed format on Amazon Web Services (AWS) servers, provided by TACC at http://aste-release1.s3-website.us-east-2.amazonaws.com/. These files are meant to be accessed with the llcreader module of the open source python package xmitgcm (Abernathey, Dussin, et al., 2020), which allows users to analyze the data without the need to actually download it. An interactive demonstration of this capability, which shows some sample calculations enabled by xgcm (Abernathey, Busecke, et al., 2020) and ECCOv4-py (github.com/ECCO-GROUP/ ECCOv4-py), is available through the Binder Project (Project Jupyter et al., 2018) at github.com/crios-ut/ aste (Smith, 2021). This repository additionally contains environment files so that any user can reproduce the computing environment necessary to analyze ASTE_R1, for instance on their own laptop.

A3. Observational constraints from ECCOv4r3 standard suite
As described in Section 2.1, the observational constraints used in ASTE include the standard ECCOv4r3 suite (Fukumori, Fenty et al., 2018) and additional high-latitude data as listed in Table 2. For a quick reference, we list the data from the ECCOV4r3 suite in Table A1 and refer the readers

Appendix B: Transport Calculation With Referenced θ/S
Here, we describe heat and freshwater transport calculations used in ASTE_R1 with respect to reference values of potential temperature (θ r ) and salinity (S r ), respectively. This serves to (1) provide calculation details for comparison to those used by previously published estimates (supplementing results presented in section 4), and 2 expose where calculation differences may prevent meaningful comparisons (following discussion in section 2.3). For budget calculations, we refer the readers to detailed descriptions provided in Piecuch (2017) and Forget, Campin et al. (2015).
In the literature, transports are often computed with nonzero referenced values θ r /S r . In section 4, we provided online transport estimates for ASTE_R1 made using non-zero references (e.g., for the heat flux through the Bering Strait). We emphasize that all offline transport calculations made using available diagnostics from ASTE_R1 (and all standard configurations of the MITgcm) will be exact only with θ r = 0 and S r = 0, as these are the values used in all online tracer equations. To support users seeking to compute ASTE_R1 transports offline assuming nonzero references, we now examine the loss of accuracy that will be incurred. This loss of accuracy depends on the amplitude of various missing terms (e.g., bolus transports and diffusive fluxes) relative to the contributions (e.g., Eulerian advection) contained in the available diagnostics. By deriving these approximations here and comparing their magnitudes with the accurate online values across important Arctic and GIN Seas gateways, we aim to identify which transports reported in Figures 16 and 17 are reliable and which ones require caution for interpretation.

B1. Accurate Transport Calculations
The horizontal transports of volume, heat, and freshwater (FW) across the Arctic Mediteranean gateways are calculated by summing the total horizontal convergence in the mass, heat, and salinity budgets, respectively, (Piecuch, 2017) as follows,  where t is the time, u E , u b the (partial-cell-weighted) ocean resolved Eulerian and unresolved bolus velocities, and u i , u sn the sea ice and snow Eulerian velocities. For each gateway across which the transports are computed, n is the normal direction at each model grid point along the transport gate and L the section length along the gate. Vertical integration is between η the sea surface and −D the ocean floor depth for volume and heat transports. Constants ρ 0 = 1029, ρ i = 910 and ρ sn = 330 are the seawater, sea ice, and snow densities in kg/m 3 . C p = 3996 J°C −1 kg −1 is the specific heat capacity of sea water, θ r the reference temperature, S r = 34.8 ppt the reference salinity, and S i = 4 ppt the constant sea ice salinity used in ASTE. θ and S are the ocean potential temperature and salinity in °C and ppt, respectively; h i and h sn are the thickness of sea ice and snow in m. F θ,dif and F S,dif are the parameterized diffusive flux of potential temperature and salinity θ r = 0 and S r = 0. For both advective and diffusive contributions to freshwater transports (Equation B3.0), vertical integration is only down to the depth of the reference isohaline  S r z .
Exact closure of heat budgets (see equations in Forget, Campin et al., 2015 andPiecuch, 2017)      .2) methods. Both methods incur errors due to reliance of the monthly S for determining the depth of the reference isohaline S r z serving as the integral limit. "off" refers to Equation B3.2 which computes the transport offline using outputs of monthly averaged Eulerian velocity u E and salinity S . "adv" and "dif" are the approximated online calculations of advective ( The variability is computed after the seasonal cycle has been removed. Note that "adv" is consistently larger than "off" (and with larger variability), but it is not possible to conclude that the online calculation is superior due to imperfect treatment of S r z . Instead, we assume higher confidence in both our FW flux estimation and our FW flux comparisons where "adv" and "off" converge.

B3. Approximations Using Monthly Mean θ, S, and u E
Lastly, we note that due to disk space and I/O restrictions, it is typical for modeling studies to save and subsequently provide only monthly averaged Eulerian velocity E u and tracers  and S for offline calculations of heat/FW transports and contents (e.g., Jahn et al., 2012;Kinney et al., 2014;Q. Wang et al., 2016aQ. Wang et al., , 2016bIlicak et al., 2016;Heuzé & Årthun, 2019). In this case, the calculation for the advective terms in heat and FW transports are further approximated due to the cross-terms involving the bolus velocity Su b and θu b being excluded: As before, inaccuracies will be incurred when S r z is determined using the monthly mean S . This is the case for all results shown for FW because no diagnostics pertinent to FW, including those of S − S r or S r z , are available standard MITgcm diagnostic outputs.

B4. Interpretation of Transports: Confidence and Caution
Figures B1 and B2 show time-series of heat and FW transports for key gateways using both the most accurate online method and approximated offline method described above. The diffusion terms for both heat and FW are at least two orders of magnitude smaller than the advection terms and can be ignored almost everywhere. The exception is at the Denmark Strait and Iceland-Faroe channels where omission of the diffusive contribution to the total heat transport leads to large errors of 30% and 100%, respectively. This shows that the estimates of tracer transports across these two gates should be interpreted with caution when computed offline using only model monthly outputs of the Eulerian velocity and tracer averages.
For FW, as all methods are approximated, the largest error is likely due to not tracking the time-evolving depth of the reference isohaline S r z . Since there is no exact calculation for comparison, it is not possible to conclude which method, "adv" using the online advective term (Equation B3.0) or "off" using the monthly mean Eulerian velocity and tracers (Equation B3.2, is "more" correct in Figure B2 and Table B1. It is likely that for gates where these two methods provide almost identical estimates (e.g., Bering, Davis and Fram Straits) we have higher confidence in our estimated FW transport. Across the CAA and the GSR the FW transport calculation depends strongly on the method employed and caution should be used in confidently reporting FW fluxes and comparing between different studies. NGUYEN ET AL.

Appendix C: Watermass Definition in ASTE_R1
Suitable specification of the characteristic salinity, potential temperature and density (S, θ, and σ) defining known watermasses can differ between observations and models due to model biases, as shown in Figure 14 in the main text for water properties in the Irminger and Labrador Seas. Watermasses can be clearly identified in ASTE_R1 as large volumes with a common formation history and distinct properties from surrounding waters, consistent with their definition in the literature. However, in regions of hydrographic bias, these watermasses will not be identified-or correctly quantified-as their observed counterparts if they are tracked following the observed values too strictly. In this appendix, we summarize the choices made in determining watermass and explore the sensitivity to these choices where appropriate. Table C1 lists the watermass properties at Fram Strait (FS) and across the Greenland-Scotland Ridge (GSR) used to identify the transports reported in Figure 15 in the main text. At the FS, the mean transports can NGUYEN ET AL.  be decomposed approximately into the West Spitsbergen Current (WSC, east of 5°E, Beszczynska-Möller et al., 2012), recirculated Atlantic Water (AW) (between 3.2°W and 5°E), and the East Greenland Current (EGC, west of 3.2°W, S ≤ 34 ppt, T ≤ 1°C). At the GSR definitions of watermasses such as the surface outflow, dense outflow, modified water, and inflow AW from Østerhus et al. (2019) and Hansen and Østerhus (2000) can be problematic when strictly applied to grid-scale average quantities. For example, the densewater in the outflow through Denmark Strait (DS) is defined in Østerhus et al. (2019) as having density anomaly σ θ > 27.8, but in ASTE_R1 outflow at the lowest depths of the strait are characterized by a lower bound of σ θ ranging between 27.28 and 27.81. For this range, the corresponding southward transports are −1.6 ± 0.9 to 0.5 ± 0.3 Sv (see Figure 15, blue color text). Similarly, over the Iceland-Faroe (IF) ridge, the southward transports of densewater defined by σ θ ≥ 22.44 or σ θ ≥ 27.55 in ASTE_R1 yield a range of −0.4 to −0.3 Sv, compared to −0.4 ± 0.3 Sv of water with σ θ ≥ 27.8 in Østerhus et al. (2019). Similar considerations applied also to dense water properties at the Faroe-Shetland (FSh) ridge (σ θ ≥ 27.81 in ASTE_R1 compared to 27.8 in Østerhus et al., 2019). For the northward flow, in addition to salinity thresholds (S ≥ [34.8,35,35.25] ppt), temperature thresholds of θ ≥ (5,4,5)°C are used in ASTE_R1 to identify the warm AW across the DS, IF, and FSh channels.

C2. Heat content of upper halocline watermass
The upper halocline watermass, defined by Timmermans et al. (2018) as a layer within lower and upper salinity bounds of S l = 31.0 ppt and S u = 33.0 ppt, respectively, was identified based on subsurface in situ observations with fine vertical sampling resolution. In ASTE_R1, with vertical grid spacing of 15-20 m within the water column depths 50-160 m, average salinity in the water column changes more abruptly than in the observations. For more accurate estimation of halocline-integrated quantities one approach is to "interpolate" the salinity in the vertical to a finer grid to find the exact depths at which salinity fits within the given bounds. Though this is often done during model-data comparisons (e.g., Grabon, 2020), the interpolation introduces additional information that was not strictly solved for by the model. An alternate approach is to vary the salinity bounds to gauge the sensitivity of the heat content within this watermass to the vertical discretization in the model. As an example, Figure C1 shows a vertical section in ASTE_R1 through the Beaufort Gyre region as defined in Timmermans et al. (2018), with the watermass bounded between a temperature maximum at depths ∼50-60 m (Pacific Summer Water, PSW, S l = 31 salinity contour) and a temperature minimum at depths ∼150 m (Pacific Winter Water, PWW, S u = 33). In ASTE_R1, negligible sensitivity is found with changes to S u , but the heat content within the upper halocline in this region changes by ∼1%-2.5% per 0.1 ppt change in S l . A change in S l of ∼0.5 ppt corresponds approximately to one depth level in ASTE_R1, and the heat content change associated with this is shown in shade in Fig AW NSF-OCE-1924546. Additional funding was provided from the ECCO project through a JPL/Caltech subcontract. Computing resources were provided by the University of Texas at Austin Texas Advanced Computing Center (TACC) and NASA Advanced Supercomputing Division at the Ames Research Center. Adjoint code was generated using the TAF software tool, created and maintained by FastOpt GmbH (http://www. fastopt.com/). Author ATN thanks M.-L. Timmermans and N. Foukal for providing heat content time-series for the Arctic halocline and eastern SPNA, and W. von Appen for the 2012-2017 processed Fram Strait mooring data used for ASTE post validation. We thank three anonymous reviewers for comments that greatly improved the manuscript.