Sea surface height anomaly and geostrophic velocity from altimetry measurements over the Arctic Ocean (2011-2018)

In recent decades the decline of the Arctic sea ice has modified vertical momentum fluxes from the atmosphere to the ice and the ocean, thereby affecting the surface circulation. In the past ten years satellite altimetry has contributed to understand these changes. However, data from ice-covered regions require dedicated processing, originating inconsistency between ice-covered and open ocean regions in terms of biases, corrections and data coverage. Thus, efforts to generate consistent Arctic-wide 5 datasets are still required to enable the study of the Arctic Ocean surface circulation at basin-wide scales. Here we provide and assess a monthly gridded dataset of sea surface height anomaly and geostrophic velocity. This dataset is based on Cryosat-2 observations over ice-covered and open ocean areas of the Arctic up to 88o N for the period 2011 to 2018, interpolated using the Data-Interpolating Variational Analysis (DIVA) method. Geostrophic velocity was not available north of 82° N before this study. To examine the robustness of our results, we compare the generated fields to one independent altimetry dataset 10 and independent data of ocean bottom pressure, steric height and near-surface ocean velocity from moorings. Results from the comparison to near-surface ocean velocity show that our geostrophic velocity fields can resolve seasonal to interannual variability of boundary currents wider than about 50 km. We further discuss the seasonal cycle of sea surface height and geostrophic velocity in the context of previous literature. Large scale features emerge, i.e. Arctic-wide maximum sea surface height between October and January, with the highest amplitude over the shelves, and basin wide seasonal acceleration of 15 Arctic slope currents in winter. We suggest that this dataset can be used to study not only the large scale sea surface height and circulation but also the regionally confined boundary currents. The dataset is available in netCDF format from PANGAEA at [data currently under review].


Sea surface height
Monthly η fields were compared to an independent satellite gridded dataset over the entire Arctic. This dataset is described by Armitage et al. (2016) and will be hereafter referret to as CPOM DOT (Centre for Polar Observation and Modelling Dynamic 120 Ocean Topography, available at http://www.cpom.ucl.ac.uk/dynamic_topography). The CPOM DOT is a regional Arctic dataset spanning the years 2003-2014, derived from sea surface height observations (relying on the Envisat and Cryosat-2 satellite missions) and a geoid model (GOCO3s). Monthly fields are provided on a 0.75°×0.25°longitude-latitude grid, up to a latitude of 82°N.
We also compared η locally, to the sum of in-situ measurements of steric height (the height component due to changes in 125 density) and bottom pressure equivalent height (related to changes in water mass). These two components were derived from temperature, salinity and ocean bottom pressure data from moorings at three sites. The moorings were located in the southern Fram Strait ([78.17º N, 0º E], hereafter FS_S), at the shelf break north of Arctic Cape, the headland of Severnaya Zemlya ([82.22º N, 94.85º E], hereafter AC), and down the continental slope north of the Laptev Sea ([78.46º-81.15º N, 125.70º (Polyakov, 2016(Polyakov, , 2019Polyakov and Rembert, 2019).

Velocity
We used measurements of near-surface velocity from two mooring lines to evaluate monthly geostrophic velocity in the Fram Strait and down the continental slope of the Laptev Sea. The Fram Strait array comprises 17 moorings located along a zonal section at 78°50' N, between the longitudes 9°W and 8°E, maintained since 1997 by the AWI (moorings F1-F10 and F15/F16; Beszczynska-Möller et al. (2012)) and the Norwegian Polar Institute (NPI, moorings F11-F14 and F17;de Steur et al. (2009)). 8°7 ' W 78°50 ' N 13 (2011-2012) Sea, data were used from four (out of six) moorings deployed in a meridional transect along the 126°E meridian within the context of the NABOS-II project (moorings M1_1 to M1_4). ADCP velocity measurements were averaged in the upper 50 m.
These deployments have been carried out twice, in 2013 and 2015. Moorings were respectively recovered in 2015 and 2018, providing a record spanning 5 years (data are available from the Arctic Data Center, Polyakov (2016Polyakov ( , 2019; Polyakov and Rembert (2019)). Moorings positions and the monthly data availability are detailed in Table 2

Finite Elements Sea ice-Ocean Model (FESOM) output
The decorrelation length scale used for the interpolation of along-track sea surface height was determined based on model data.
We used monthly geostrophic velocity derived from sea surface height outputs of the Finite Elements Sea ice-Ocean Model (FESOM) version 1.4. FESOM 1.4 is a coupled sea ice -ocean model, working on a triangular unstructured mesh. The model has a resolution of 4.5 km over the Arctic Ocean and has been described and validated by Wang et al. (2018) and Wang et al. 160 (2019). The run used in this work is the historical run described by Wang et al. (2020), forced by atmospheric reanalysis data of JRA55-do v.1.3 (Tsujino et al., 2018).

Methods
In this section we provide a description of the in-situ data processing (Sect 4.1) and of the steps followed to derive altimetry monthly fields from along-track satellite measurements (Sect. 4.2,4.3 and 4.4). In the last two subsections we provide details 165 on how we performed the comparisons with independent datasets and how we computed the seasonal cycle (Sect. 4.5, 4.6).

Processing of in situ data
In-situ steric height anomaly (η S ) and bottom pressure equivalent height anomaly (η P ) were computed from measurements of water density and ocean bottom pressure. The relation between η and the time anomaly of i) the vertical density profile (ρ (z)) and ii) the ocean bottom pressure (P b ), is derived by integration of the hydrostatic balance from the sea surface down to the 170 bottom depth, D (Eq. 5): where g is the gravitational acceleration and ρ 0 is a reference ocean water density, set to 1028 kg m −3 . Based on this relation, we defined η S and η P as (Eq. 6): We computed η S and η P from measurements of ρ (z) and P b , at the FS_S, AC and M1_4p6 mooring positions. Vertical profiles, ρ (z), were obtained from temperature and salinity profiles using the UNESCO (1983) formula for density. In turn, temperature and salinity profiles were obtained from moored-sensor data by linear interpolation on a regular pressure grid (2 dbar) between the shallowest sensor (FS_S = 50 m, AC = 50 m, M1_4p6 = 26 m) and the deepest sensor (FS_S = 729 m, AC = 1448 m, M1_4p6 = 700 m). Data was extrapolated from the shallowest sensor to the sea surface assuming constant temperature 180 and salinity, and equal to the uppermost measurement. Below the deepest sensor we assumed the density anomaly to be zero and did not perform extrapolation to the bottom. This conservative approach might have resulted in the underestimation of η S .
Ocean bottom pressure records P b were de-tided using the Matlab function t_tide (Pawlowicz et al., 2002) and instrumental drifts were removed. Unfortunately the time series at FS_S exhibited large pressure anomalies, developing on timescales of several months, whose amplitude was at least one order of magnitude too large to be explained by changes in ocean currents.

185
Therefore, we high-pass filtered this time series with a cutoff frequency of 2 months. All other bottom pressure time series were not affected.

Along-track sea surface height anomaly
We generated an Arctic-wide dataset of along-track η by merging the AWI and RADS η datasets. Inconsistencies between the two datasets were reduced by: i) creating a uniform along-track sampling, ii) reducing biases due to different retracking 190 algorithms, and iii) substituting geophysical corrections where two different corrections were used in the two source products.
Here we first give details about these methods and present an estimate of the η observational uncertainty at the end of the section.

Merging leads and open ocean data
Prior to merging the AWI and RADS datasets we standardized their along-track sampling rates, which originally were respec-195 tively 300 m and 7 km. With this aim, the AWI dataset was first smoothed with a moving window of 7 km and then linearly interpolated, following time, onto equally spaced locations (7 km) along the satellite tracks.
A step-like variation in the η observations at ocean-ice transitions appeared because different models are used to retrack signal returns in ice-covered and ice-free regions (Fig. 2). This is commonly referred to as the "lead-open ocean bias" . Due to the technical nature of this bias, it is difficult to determine the true bias in this post processing phase.

200
Therefore, we simply removed the offset in η between leads and open ocean data, estimated directly from the along-track η .
To do so, we identified along-track transitions between the AWI and RADS datasets where the gap was shorter than 200 km.
Then we selected observations from the two datasets within an along-track distance of 200 km from the last ice-covered data   (Fig. 2a). We calculated the offset at each transition as the difference between the mean values of the selected ice-covered and ice-free η observations. The offset distribution is different by season (Fig. 2b), but no sensitivity was shown to the ocean-205 ice transition direction or to the gap between last ocean and first ice data points. Accordingly, we derived monthly offset values as the median of the offsets in that month.

Corrections
As second step, we checked that all corrections applied to the satellite range R (Eq. 1) were consistent between ice-covered and ice-free regions (Table 3 lists the products used here). Standard corrections (European Space Agency, 2016) were applied to 210 both regions to account for i) the reduction in satellite signal speed caused by the presence of the atmosphere (dry gases, water vapour, ions); ii) the difference in reflection properties of wave troughs and crests at the sea surface (sea state bias correction, applied solely in the open ocean); and iii) solid earth tides and ocean tides.
A further correction removes the high frequency ocean response to atmospheric pressure and wind forcing. For ice-covered regions ESA suggests using an Inverted Barometer (IB) formula. Here instead, we applied the Dynamic Atmosphere Correction 215 (DAC, Carrère and Lyard (2003); Carrère et al. (2016)) to both, ice-covered and ice-free regions. DAC is conventionally used in the global ocean because it better suppresses high frequency variability, avoiding aliasing of sub-monthly temporal changes into spatial variability (Carrère and Lyard, 2003;Quinn and Ponte, 2012;Carrère et al., 2016).
No study to date shows which of the DAC and IB corrections performs better in ice-covered regions. However, gaining insight into this issue was relevant to us because residual sub-monthly variability emerges in the monthly η fields as meridionally 220 elongated patterns (meridional "trackiness", Stammer et al. (2000)). Therefore, to support our choice of using DAC over IB,   (2005) Ocean Tide FES2004 Lyard et al. (2006) Solid Earth Tide Cartwright model Cartwright and Edden (1973) Geocentric Polar Tide Instantaneous Polar Location files (sourced from CNES) Wahr (1985) Orbit GDR-E European Space Agency (2016) we looked at which of them reduced the η standard deviation the most with respect to the uncorrected η (see Appendix A).
Results showed that DAC outperforms the IB in shallow shelf regions (particularly the East Siberian Sea and the Chukchi Sea) and that they perform equally well over the deep basins ( Fig. A1). For instance, in the East Siberian Sea the DAC reduced the uncorrected η standard deviation by 50% at periods shorter than 20 days, in contrast to no reduction when applying a simple 225 IB (see Table A1). The improvement of DAC with respect to IB over the shelves appears also in the η monthly grids, where meridionally oriented patterns of η are evidently reduced (two examples are given for the months of November 2014 and November 2017 in Fig. A2).

Final along-track dataset and uncertainty estimate
Following the steps above we generated a final merged along-track dataset, composed of two sub-datasets; one for the ice-230 covered region and one for the ice-free region. The consistency of these two sub-datasets is indicated by their comparable Arctic-wide standard deviation over the period 2011-2018 (10.4 cm and 10.5 cm, respectively). An example is illustrated for the month of July 2015 in Fig. 3, which shows a smooth transition between the two sub-datasets. We note though, there is some residual sub-monthly variability. For instance, Fig. 3a shows a decrease of η of ∼20 cm north of Greenland between the first and the fourth week of July 2015. The residual sub-monthly variability is one of the two main contributions to the error on the monthly η fields. This contribution is estimated as part of the processing algorithm during the interpolation phase (Sect.

4.3.3).
The second contributor to the error is the observational uncertainty.
The observational uncertainty associated to the along-track η stems from several sources, namely the altimeter measurement uncertainty, the waveform retracking method, the corrections and orbit uncertainty. Given the difficulty of assessing the contribution of each of these sources, we provide here a comprehensive estimate of the observational uncertainty based on the 240 differences of the along-track η at satellite tracks crossovers (Fig. 4). We first defined crossovers as those pairs of η observations within a distance of 7 km. We excluded pairs belonging to the same satellite pass by verifying that they are separated by more than one hour. Considering the large number of data, we organised observations in an equal area grid of about 100 km  and computed the η differences at crossovers only within selected cells (red dots in the inset of Fig. 4). We expect that the η differences tend to the observational uncertainty as the crossovers get closer in time. Thus our uncertainty estimate was given 245 by the average η difference at crossovers separated by no more than 5 days, which is 4.2 cm.
This analysis provides additional information about the η de-correlation time scale. The η crossover difference increases with time above the uncertainty due to local variability. Fig. 4 shows that variability increases by ∼2 cm in the first 20 days, then by a further ∼1 cm after 6 months, until it finally reaches a plateau. This indicates that, at time scales shorter than one year, η has a short de-correlation time scale of 20 days and a long de-correlation time scale of 6 months.

Gridded sea surface height anomaly
We produced monthly η fields over the period 2011-2018, by interpolation of the along-track data onto a longitude-latitude grid of resolution 0.75°×0.25°, from 60°N to 88°N. Gridding creates regular fields from irregularly distributed data points.
Below we first describe the technique used to reduce residual sub-monthly variability in the along-track η (Sect. 4.3.1). Then give details about the interpolation method used, including the selection of the length scale and signal to noise ratio (Sect.

Minimisation of sub-monthly variability
The residual sub-monthly variability (described in Sec. 4.2.3) produces marked meridional trackiness if the interpolation is performed on a monthly set of η observations (see Fig. 5a, July 2015). To reduce this variability we performed the interpolation on weekly data subsets instead. The monthly maps were then obtained as the average of four weekly maps, with the associated 260 error given by the quadratic sum of the weekly error maps (shown below in Fig. 6). Fig. 5b shows data along a latitude circle as an example of the trackiness reduction obtained thanks to this approach.

Interpolation via DIVA
In oceanography, regular fields are commonly created via a technique called optimal interpolation (Bretherton et al., 1976). This approach though, has the disadvantage of generating the interpolated field based on an a priori functional form of the covariance between data points. The variational inverse method that we use here is an alternative which does not assume an a priori behaviour of the analysed field. Instead, it improves the quality of the interpolated field (e.g., closeness to data, smoothness), by minimizing a cost function dependent on the data and few input parameters (Troupin et al., 2012). An additional advantage of this method is that it decouples ocean basins separated by land by minimizing the cost function on a finite elements mesh, 270 whose nodes are not connected across land.
We used a specific realization of the variational inverse method, namely the Data-Interpolating Variational Analysis (DIVA, Troupin et al. (2012)). Further information is provided in Appendix B. As input to DIVA, the along-track η data and two parameters, the length scale L and the data signal to noise ratio λ are provided. DIVA gives the possibility to introduce an anisotropic weighting of the data points, which we tested via an input vector field.

275
The selection of values for λ and L was done as follows. First we carried out twenty interpolation runs with λ =0.1, 0.3, 1, 10 and L =20 km, 50 km, 100 km, 300 km, 600 km. These values of λ cover the cases when a) data contain little information with respect to the noise level or the data are not representative of the monthly mean (λ =0.1), and b) when the noise level is only 10% of a real signal (λ =10). The range of length scales L covers Arctic Ocean circulation scales, from mesoscale dynamics to basin wide flow (e.g., Nurser and Bacon, 2014;Armitage et al., 2016). Since λ and L are the main input parameters for 280 the interpolation, we did not include the anisotropic weighting initially, but its effect was tested separately at a later stage (described below).
Next, we compared the results of each interpolation run with two independent datasets: η P + η S time series at the FS_S, AC and M1_4p6 moorings (Sect. 3.2.1); average geostrophic transport from the FESOM model outputs (Sect. 3.3), normal to two transects; one crossing the 285 opening of the Laptev Sea into the Eurasian Basin and another across the north-western Fram Strait (both transects are indicated in Fig. 1).
The aim of these comparisons was to select the λ and L values whose interpolated field best captured i) the high temporal resolution patterns exhibited by the in-situ data and ii) the large scale features simulated by the model.
The comparisons were assessed by defining a score S(λ, L), such that (Eq. 7): where C i (λ, L) represents the correlation of a given run (with parameters λ and L) with the in-situ or model time series i; RMSD i (λ, L) is the root mean square difference between each run and the time series i; and X i,max and X i,min are the maximum and minimum values of the interpolation (either η or (u g , v g )) in each comparison. This score is designed to be  weight.
The scores (S) are presented in Table 4. They vary only by a small fraction, indicating that the patterns in the altimeter fields are insensitive to the choice of parameter values. There is though, a minimum score for λ = 1 and L = 50 km, which were thus selected as a suitable set of parameter values.
Finally, we tested the anisotropic weighting (Eq. B3) by including it in interpolation runs with the selected λ and L. As vector 300 field we used the long-term mean (1992-2012) geostrophic velocity field from the FESOM model. We used two amplification factors a (a = 100, a = 10 6 ) to represent the effect of weak or strong advection. Results showed the effect of advection of η by the mean flow to be negligible, even with high a. Therefore, we did not apply this constraint to compute the final η gridded fields.
The error maps associated with the interpolation were provided by DIVA through the poorman's estimate method (Troupin 305 et al., 2012). This method circumvents the high computational cost of calculating the real covariance in DIVA by assuming a constant covariance between data points. The poorman's estimate method generates maps of relative error, given as fraction of the variance of the background field (Troupin et al., 2012). These maps allow to assess the data coverage given by the distribution of the data in space scaled by the decorrelation scale L. An approximate way to scale the relative error to observational units is by multiplication with the standard deviation of the input data. For the final product we provide the relative error maps.

Error on monthly fields
The error in the monthly η fields comprises a component arising from the observational uncertainty and another arising from the sub-monthly variability. The error arising from the observational uncertainty varies spatially depending on the data distribution and the interpolation method. However, given that the poorman's estimate method does not provide an estimate of the absolute error, we used an alternative procedure to yield an average estimate of the standard error components over

Gridded geostrophic velocity
The geostrophic velocity was computed on the output grid following Eq. 4, with partial derivatives approximated by finite differences. The components of velocity on the longitude-latitude grid at indices i, j are given by (Eq. 8): where θ and Φ are latitude and longitude converted to radian angles, f = 2Ωsin(θ) is the Coriolis parameter (with Ω = 330 7.29 · 10 −5 s −1 ) and R e is the Earth radius. In the equation above, η i,j is the η field at indices i, j obtained by adding the gridded η to the DTU17MDT linearly interpolated to the grid.

Comparison to independent datasets
The monthly η maps were first evaluated over the whole Arctic Ocean by comparison to the CPOM DOT product. The comparison was done at grid points south of 82º N (the northernmost latitude covered by the CPOM DOT), for the period 335 January 2011 to December 2014. Both datasets were referred to their own temporal average over this period.
In the next step we evaluated our η fields locally via comparison with time series of in-situ sea surface height η i = η P + η S (Eq. 5 and 6) from moorings FS_S, AC and M1_4p6. The altimetry η was linearly interpolated to the in-situ locations, and the time average over the period of mooring deployments were removed from in-situ and altimetry observations.
Lastly, we assessed geostrophic current fields (u g , v g ) by comparison to in-situ measured currents from the two mooring lines 340 crossing the Fram Strait and the Laptev Sea continental slope. We compared the (u g , v g ) component normal to the transects, linearly interpolated to the moorings positions, to monthly averages of the in-situ measured velocities normal to the transects.
These are hereafter referred to as v n and v ni respectively, positive northward in the Fram Strait and eastward in the Laptev Sea.
The comparison was limited to those mooring locations where more than 24 months of in-situ data were available at the time of manuscript preparation.

345
To establish spatial scales over which altimetry-derived currents approximate best the in-situ measured currents, we compared v n and v ni spatially averaged over different sets of neighbouring moorings. The averaging was performed by assigning to each mooring a weight proportional to its distance to the two neighbouring moorings (e.g., for mooring j: w j = dj−1,j +dj,j+1 2 , where d is the distance). We performed three tests in the Fram Strait and two at the Laptev Sea continental slope, using moor- we evaluated the percentage of variance explained by the seasonal to interannual frequency band (lower than 4 months) and by the intra-annual frequency band (higher than 4 months) as (Eq. 9): where x is the horizontally averaged velocity time series (tests 1 to 5), and x F is the correspondent filtered time series.

Seasonal cycle
The seasonality of the Arctic sea level and surface currents has been studied in several previous works (e.g., Volkov et al., 2013;Armitage et al., 2016;Beszczynska-Möller et al., 2012;Baumann et al., 2018), giving us the opportunity to assess our dataset based on this literature. We defined the seasonal cycle of η , following Volkov et al. (2013), as the harmonic least-square fit to η with period of one year (Eq. 10): where t is the number of the month in the time series (t = 1 correspond to January 2011) and P = 12 is the oscillation period.
We evaluated the fraction of variance explained by η seas at each grid point following Eq. 9, with η as x and η seas as x F .

Results
Here we first describe the characteristics of the η and geostrophic velocity (u g , v g ) monthly maps, then show the results of 365 their comparison against independent datasets, and lastly present the η and (u g , v g ) seasonal cycle.

Monthly fields of sea surface height anomaly and geostrophic velocity
Given our data set spans 96 months within the 2011-2018 period, here we present results from the month of July 2015 as an example to describe general characteristics of a given map. Fig. 6 shows fields of η , relative error (associated with the interpolation) and (u g , v g ) for July 2015. The description below makes reference to the Arctic Ocean sub-regions and long 370 term mean surface circulation pathways presented in Fig. 1.
In the η monthly fields we find that there are extended regions of either positive or negative values. From Fig. 6a for instance, it can be appreciated that η is positive in the Nordic Seas and across the Arctic Deep Basins, but negative over the Eurasian Shelves. η also varies within these regions, being maximum (∼10 cm) north of 85º N, and minimum in the Laptev Sea. Superimposed on these large scale patterns, residual meridional trackiness appears south of 80°N, especially in shallow 375 areas, where the error related to sub-monthly variability is relatively high (Fig. 5c). For example, enhanced trackiness is visible in the Barents Sea, where previous work has shown that intra-seasonal sea surface height variability explains between 50% and 80% of the total variability (Volkov et al., 2013).  (Fig. 6b). Largest relative errors are found in regions with data gaps (see data 380 distribution in Fig. 3a): i) south of 75°N, where the distance between the satellite tracks increases considerably; ii) in a zonal band around 80°N, where the weekly data distribution is not uniform due to the satellite orbit geometry; and iii) in the Canada Basin, where there is a large and persistent marginal ice zone. Neither the AWI nor by the RADS processing provide data over marginal ice zones.
In Fig. 6c we present the geostrophic vector field (u g , v, g ), with background colors highlighting monthly speed anomalies Greenland Shelf). The root-mean-square deviation (RMSD) is presented in Fig. 7b, showing low values (2 cm to 4 cm) over 75% of the domain, including most of the regions with water depth greater than 100 m. The RMSD is high (7-8 cm) over the East Siberian Sea and Chukchi Sea, where the sub-monthly variability is most enhanced.
Time series of the comparison of η with in-situ data are shown in Fig. 8, and the associated RMSD and correlation coefficients are presented in Table 5. The correlation between the η and η i time series is relatively low (0.3 to 0.5), but significant 405 (p-value < 0.05). η and η i follow roughly a similar pattern, varying within a range of ±10 cm over the comparison period at all three sites. At the FS_S mooring there are hints of a seasonal oscillation, with the signal decreasing from October 2016 to March-April 2017 and then increasing towards October 2017 (Fig. 8). At the AC and M1_4p6 moorings short term variability appears in phase at times, for instance between December 2016 and May 2017 in the former, and between September 2014 and February 2015 in the latter (Fig. 8). There are however, differences as large as the altimetry variability in some months, which 410 is reflected in the RMSD and the η standard deviation (Table 5).  Table 5. Comparison between altimetry and in-situ sea surface height anomaly. The first row show the correlation (Pearson's correlation coefficient, p-values were computed using the effective number of degrees of freedom (Emery and Thomson, 2001)). The second row shows the RMSD between altimetry η and in-situ η i at moorings FS_S, AC and M1_4p6.

Velocity
The agreement of in-situ and altimetry-derived velocities is summarized in Table 6, which presents their correlation and RMSD, together with the mean velocities and standard deviation at each mooring (computed over the months when in-situ data were available). Time series across the Fram Strait and Laptev Sea mooring lines are displayed in Fig. 9 and 10. In the Fram Strait, 415 the correlation is significant (p-value < 0.05) and higher than 0.3 at moorings F2 to F5. At these 4 moorings, both the mean v n and v ni are consistently positive and comparable or higher than the corresponding standard deviation. The correlation was highest at mooring F3; the mooring with the longest continuous time series. In the Laptev Sea continental slope the correlation is highest at the M1_1 mooring. Further down the slope the correlation is lower, being still significant at mooring M1_3, but non-significant at moorings M1_2 and M1_4.

420
There are differences between altimetry and in-situ data in terms of spatial and temporal resolution. The mean v n shows low spatial variability and smooth transitions between nearby sites ( Table 6). Note that this variability is governed by the averaging scales underlying the DTU17MDT product. The scales captured by the DTU17MDT are defined by the resolution of the geoid model used to compute it. Previous studies indicate that these scales are not smaller than 100-150 km (Gruber and Willberg, 2019;Bruinsma et al., 2014;Farrell et al., 2012). These large scales contrast with the high spatial variability of the v ni mean 425 flow, which is derived by pointwise measurements. This is shown for instance by abrupt changes between moorings F15 and F7 (50 km apart) or between M1_1 and M1_2 (11 km apart). High spatial variability observed by the mooring data is ascribable to the small Arctic first baroclinic Rossby radius, which is below 10 km in the two study regions (Nurser and Bacon, 2014; Pnyushkov et al., 2015). Furthermore, the time variability of the mesoscale processes is smoothed out in the altimetry dataset due to the 50 km decorrelation scale applied through the interpolation. This is reflected in the v n standard 430 deviation, which is about four to five times smaller than that of v ni at most moorings.
Tests 1 to 5 show the scales over which spatial averaging improved the comparison (Table 7). In the Fram Strait, averaging over moorings F3 to F5 (test 2) yielded a correlation higher than that using data only from the F3 mooring (Tables 6 and 7).
Results from tests 1 and 3 yielded correlations comparable to that from F3. All three tests reduce the RMSD by about 2-3 cm with respect to that at F3. At the Laptev Sea continental slope, neither test 4 nor test 5 improved the correlation with respect to 435 the comparison at the M1_1 mooring. Both tests though reduce the RMSD with respect to the one at M1_1 (2-4 cm lower).
With regards to temporally filtering of the time series in tests 1 to 5, we find that seasonal to interannual frequencies explain most of the variability both in v n and in v ni . They constitute about 80% of the total variability in the Fram Strait, and about 90% at the Laptev Sea continental slope. In this frequency band v n and in v ni correlate better or equally than without filtering (Table 7), whereas in the intra-seasonal frequency band the correlation worsen.

440
The fact that seasonal to interannual frequencies explain a high percentage of the total variability can be attributed to the dominant seasonal oscillations. These are visible for instance in Fig. 9a and Fig. 10a

Seasonal cycle
In the following text we give an overview regarding the seasonal cycle observed in η and (u g , v g ), highlighting which are the regions where it explains a high fraction of the total variability.

Sea surface height anomaly
The amplitude A and the phase α of the η seasonal cycle (i.e., month when the maximum occurs) are shown in Fig. 11. Regions where the seasonal cycle explains less than 20% of the total variability are blanked in Fig. 11b. The seasonal cycle explains more than 20% of the total variability in shallow shelf regions (with a peak of 60% in the Barents Sea) and in few deep regions like the southwestern Canada Basin, the Baffin Bay and the Nordic Seas. Here, the amplitude ranges between 3 cm and 8 cm 455 (Fig. 11a). Instead, seasonal variability seems to play little role in the Arctic Deep Basins, where its amplitude is <1 cm. η seas is maximum in early winter across the Arctic Ocean, even though not uniformly. η seas peaks earliest (September-October) in the Nordic Seas, the Chukchi Shelf and the Baffin Bay (Fig. 11b). A maximum around November-December is found further inside the Arctic, in the Barents and western Kara Seas, in the Laptev Sea and East Siberian Sea and in the southwestern Canada Basin. The maxima in the eastern Kara Sea (January) and the northeast Greenland Shelf (January-February) are later 460 in the winter.
In Fig. 11c we also display the observed η monthly climatology in selected regions, computed as the January to December monthly averages over the years 2011-2018. We see that the harmonic fit is a good approximation of the climatology in most of these regions, with few exceptions. For instance, in the Canada Basin, the Baffin Bay and the northeastern Greenland Shelf the climatology exhibits a secondary peak in April.  Table 8. From the comparison between 470 summer and winter current anomalies we observe a basin wide, coherent seasonal acceleration of the Arctic slope currents in winter and a deceleration in summer. The speed of these slope currents peaks between September and April. Namely, currents along the Nansen Basin shelf break, between the Fram Strait and the Lomonosov Ridge, peak in early winter (September to December); currents along the eastern shelf break of the Nordic Seas, in the Barents Sea and in the Baffin Bay peak in mid winter (November to February); the East Greenland Current peaks in late winter (February to April). Seasonality is also 475 recognisable in some currents not along the continental slopes. For instance, currents along the Siberian coasts, in the Kara Sea

480
The dataset presented in this paper provides 8 years of monthly maps of altimetry-derived sea surface height anomaly η up to 88°N. In addition, we also provide the associated geostrophic velocity (u g , v g ), which was not available before north of 82°N . We performed an Arctic wide comparison against an independent altimetry product. Results from this comparison indicate isolated areas where the correlation between data sets is low. Thus, here we discuss whether this is related to the methods used. We also carried out a comparison to in-situ data of both sea surface height and surface velocity in three seasonally ice-485 covered regions. We discuss these results in terms of the spatial and temporal resolution of our altimetry-derived velocity and the underlying dynamic regimes. Finally we put our findings on the seasonal cycle of sea surface height and geostrophic flow in the context of previous literature.

Impact of methodology
The comparison with the CPOM DOT (Sect. 5.2.1) yielded a correlation higher than 0.7 over most of the domain, but lower in 490 some regions, with non-negligible differences between the datasets there. What are the methodological steps that may generate differences between these two data sets?
Regional differences might occur due to different data density, which results from distinct algorithms used for the processing of waveforms. In our comparison the correlation is low in some areas of the ice-covered Arctic. In ice-covered regions the detection of leads is based on surface classification techniques. These differ substantially between studies, and are to date a 495 source of uncertainty (Dettmering et al., 2018). For instance, more observations are discarded the more conservative a technique is. This yields lower uncertainty, but also lower data density in the central and western Arctic, where the most compact multi year ice is located and lead density is low (Willmes and Heinemann, 2016). Thus, differences in the ice-covered regions are expected, given that Cryosat-2 observations are not classified in the AWI and CPOM DOT datasets using the same parameters and thresholds (Ricker et al., 2014;Armitage et al., 2016).

500
Generating data over the marginal ice zone still represents a challenge to overcome. This is because neither ocean-type retrackers nor ice-type retrackers are suitable to process altimetry waveforms there. Thus the coverage of these areas depends on the integration of data from ice-free and ice-covered areas, e.g. via interpolation, which is thus less constrained by actual data. It is perhaps not surprising then, that the correlation lower than 0.7 resulting from our comparison corresponds to open ocean areas of the central Arctic, where large patches of low ice concentration form at the end of summer. While we are aware 505 that in our case neither the AWI nor the RADS dataset provide data in the marginal ice zone, not enough information on the CPOM DOT data coverage is available from Armitage et al. (2016).
We also found that sub-monthly η variability in the Arctic can yield substantial noise in the monthly maps, especially on the shelves (Sect. 4.3.1). To reduce this noise we took two steps: first, we substituted the IB with the DAC correction in ice- This approach yielded improvements relative the CPOM DOT dataset. For example, in the East Siberian Sea and the Chukchi Sea regions, where the sub-monthly variability is strongest, our monthly η fields have average spatial standard deviation of 6 cm for the period 2011-2014. This significantly lower than the corresponding value of 11 cm for the same region and period in the CPOM DOT. This might be attributable to larger unresolved sub-monthly variability in the CPOM DOT.
6.2 Comparison between satellite altimetry retrievals and in-situ data 515 Independent in-situ sea surface height from mooring data were used to assess our altimetry product in two separate regions of the central Arctic, i.e. the Fram Strait and the Nansen Basin. Results showed that altimetry and in situ data yield roughly consistent patterns, e.g. a clear seasonal signal in the Fram Strait and enhanced monthly variability in the Nansen Basin. Thus, both datasets consistently suggest that the sea surface height variability has different dominant time scales in the two basins.
Correlations for our open ocean comparison between altimetry and mooring observations were low relative to previous studies 520 which compared altimetry and near-shore tide gauge measurements (Volkov and Pujol, 2012;Armitage et al., 2016;Rose et al., 2019). This can be expected given that tide gauges measure sea surface height, directly comparable to altimetry. Also, we expect sea surface height variability near the coast to show larger amplitudes than in the open ocean (see Fig. 11 consistent with other studies comparing altimetry to in-situ observations. RMSD values in the range of 2 to 12 cm has been found both from the comparison with tide gauges across the Arctic (Volkov and Pujol, 2012;Armitage et al., 2016;Rose et al., 2019) and with steric height measured in the Arctic Deep Basins (Kwok and Morison, 2011).
One comparison of altimetry-derived currents with moored currents velocity was done previous to this work by Armitage et al. (2017). Despite this comparison was performed in a region other than those considered here, results are consistent with 530 those of this work. Their correlation to ADCP measurements in the interior of the Beaufort Sea, lower or equal to 0.54, are in line with our findings at most mooring sites except for mooring M1_1, which shows higher correlation (0.77). RMSD of 1-2 cm s −1 over currents of 1-6 cm s −1 also agrees well with the RMSD that we find in the interior of the Eurasian Basin (1.3-1.9 cm s −1 ).

Temporal and spatial resolution of altimetry-derived currents 535
The large spatial extent of mooring velocity measurements and their long deployment period allowed us to examine the correlation of altimetry and in-situ velocity over both different dynamic regimes and spatio-temporal scales.
We found that correlation was higher in regions where the flow variability is dominated by steady currents (i.e. boundary currents) and lower where it is dominated by nonstationary eddy activity. In the Fram Strait, altimetry and in-situ data had the highest correlation on the shore and continental slope east of 5°E, within the West Spitsbergen Current, with maximum is strongest (Aksenov et al., 2011;Baumann et al., 2018). On the contrary, in both regions the correlation broke down where mean currents are slow and the mesoscale activity is enhanced. Namely, correlation was low and non significant at moorings in the central Fram Strait, where the surface circulation is dominated by westward eddies propagation 545 Hattermann et al., 2016). Similarly, correlation was low in the offshore part of the Laptev Sea continental slope, where current speed is low and eddy activity increases (Pnyushkov et al., 2015Baumann et al., 2018). The correlation varies with the dynamic regime due to the different sampling of mesoscale activity by moorings and by altimetry. Mesoscale features are not resolved in our monthly altimetry fields because of the 50 km smoothing scale used in the interpolation. This is equivalent to about ten times the local first baroclinic Rossby radius (Nurser and Bacon, 2014;von Appen et al., 2016;Pnyushkov et al., 550 2015), which roughly sets the horizontal scale of mesoscale eddies.
In-situ surface velocities were used thus to evaluate the effective spatial and temporal resolution of altimetry-derived currents.
In the region of the West Spitsbergen Current, the correlation was higher when averaging horizontally the in-situ observations over 50 km relative to 20 km. This indicates that the boundary current variability as observed by our altimetry-derived velocity agree most closely with the in-situ observed variability when both are averaged across at least 50 km. Slightly lower correlation 555 was shown when averaging data further into the central Fram Strait (see test 3 in Table 7), where we entered a different dynamic regime. Eddies are there a source of variability at intra-seasonal time scales, which is not resolved by altimetry and which biases the large scale average velocity from moorings. The considerations above suggest that our monthly geostrophic velocities can resolve seasonal to interannual variability of boundary currents wider than about 50 km. Mesoscale intra-seasonal variability is instead not resolved.

Seasonality
The sea surface height seasonal cycle is driven by changes in its steric component (due to sea ice melting and refreezing, solar insolation) and mass component (due to water accumulation or release, precipitation, evaporation, river runoff). Previous studies identified the seasonal cycle as the dominant component of the sea surface height variability in the Arctic (e.g., Volkov et al., 2013;Armitage et al., 2016;Müller et al., 2019). Our results confirm these findings, showing that this variability explains 565 a fraction higher than 20% of the total variability in large areas of the Arctic, including the Arctic Shelves, the Nordic Seas and part of the Canada Basin. Additionally we found that seasonal to interannual variability explains more than 80% of the geostrophic velocity variability within boundary currents in the Fram Strait and at the Laptev Sea continental slope.
Large scale features emerge in the seasonal cycle of η and (u g , v g ). For instance, η has over most of the Arctic seasonal maximum in the winter months between October and January. This is in agreement with previous studies of both steric height 570 and ocean mass seasonality from in-situ data. From hydrographic profiles the steric height was found to peak between September and November in the Greenland and Norwegian Seas (Siegismund et al., 2007), in the central Barents Sea (Volkov et al., 2013) and in the Canada Basin (Proshutinsky et al., 2009) documented for several regions by previous studies based on in-situ data, satellite data and model output (Table 8).
Our dataset is thus able to describe the seasonality of sea surface height and geostrophic currents across the Arctic consistently with previous studies.

Conclusions
With this work we aim to support and contribute to basin scale observational studies of the Arctic Ocean circulation by 585 providing a new gridded product of satellite-derived sea surface height anomaly (η ) and geostrophic velocity (u g , v g ). In this paper we present Arctic-wide monthly maps of η and (u g , v g ), spanning the years 2011 to 2018, covering both the ice-free and ice-covered parts of the ocean. We believe that this well documented and validated Arctic-wide dataset will be of help to the scientific community to further understand the Arctic Ocean surface circulation and sea surface height variability down to 50 km scale at seasonal to interannual time scales.

590
In our description and discussion of processing methods we find that residual sub-monthly variability in the Arctic Ocean is a source of noise for the η monthly maps. Therefore, we average four weekly maps of interpolated data. Further, the integration of altimetry data over ice-free and ice-covered regions raises limitations in terms of data coverage between these two regions.
Isolated differences in the comparison with the CPOM DOT are also attributable to the correction of sub-monthly variability and the data coverage in ice-covered regions, even though there is overall good agreement.

595
The comparison to in-situ sea surface height and near-surface velocity shows that the agreement varies between regions depending on the nature and scales of the variability. Geostrophic currents and in situ velocity have the highest correlation in regions where a stable flow (e.g. boundary currents) dominates the mesoscale eddy activity. There, the correlation is improved by spatially averaging in-situ data over cross-flow distances of at least 50 km and by filtering out intra-seasonal variability.
Overall, the comparison with in-situ data yields correlation and RMSD consistent with previous studies. The average cor-600 relation of η with the sum of steric height and bottom pressure equivalent height is 0.41 and the average RMSD 4.1 cm. The correlation between (u g , v g ) and near-surface moored velocity is highest at mooring sites within boundary currents both in the Fram Strait (0.53) and at the Laptev Sea continental slope (0.77). The average RMSD between velocities is 6.5 cm s −1 in the Fram Strait and 3.4 cm s −1 at the Laptev Sea continental slope.
Large scale patterns emerge from our preliminary analysis of the seasonal cycle of η and (u g , v g ). The η shows a basin 605 wide, coherent seasonal cycle, with maximum between October and January, and higher amplitude on the shelves. The (u g , v g ) features intensification of the Arctic slope currents in winter and weakening in summer. These characteristics are in agreement with several regional studies of the Arctic Ocean sea surface height and boundary currents. The DAC corrects the local and the dynamic ocean response (waves) to pressure and wind changes and is derived from the sea surface height output of a barotropic model (Carrère and Lyard, 2003;Carrère et al., 2016). Up until the early 2000s, the effect of atmospheric pressure and winds on sea surface height had instead been corrected using an Inverse Barometer formula (IB, e.g., Ponte and Gaspar (1999); Carrère and Lyard (2003)). In the IB assumption, the sea surface height responds locally to changes in pressure, decreasing of approximately 1 cm for each increase in pressure of 1 mbar (atmospheric loading). However,

630
it has been shown that the IB is not always a good approximation of the ocean response, especially on time scales shorter than 20 days (Carrère and Lyard, 2003).
Studies from the last two decades have shown that the deviation of ocean barotropic response from a simple IB is larger at higher latitudes (Stammer et al., 2000;Vinogradova et al., 2007;Quinn and Ponte, 2012). For instance, Quinn and Ponte (2012) found that the coherence between ocean mass variability (from GRACE satellite) and altimetric sea surface height variability 635 after applying only IB correction, increases with latitude. In the Arctic, the effect of pressure and wind forcing is not only local but also travels eastwards over the shelves in the form of mass waves (Fukumori et al., 1998;Danielson et al., 2020;Fukumori et al., 2015;Peralta-Ferriz et al., 2011). However, to date there is no study showing the effect of this waves on sea surface height measured from altimetry.
Here we compare the reduction in altimetry standard deviation obtained by applying DAC with respect to IB in ice-covered To understand which frequency bands have mostly contributed to this improvement, we take as an example the East Siberian Sea (yellow square indicated in Fig. A1a). We generated three time series of uncorrected η , η corrected by IB and η corrected 645 by DAC, averaged with timestep of 1 day over the indicated region. For each year we analysed periods between November and July, which are the only months when data from leads are available. For each time series, we computed the standard deviation in frequency bands with periods T > 20 days, 5 days < T < 20 days, T < 5 days (Table A1). Results show that DAC reduced the uncorrected η standard deviation by 50% at periods shorter than 20 days, in contrast to no reduction when applying a simple IB.

650
Furthermore, standard deviation at periods between 20 days and 5 days is larger then 60% the standard deviation at periods longer than 20 days, confirming how high frequency variability is relevant in the Arctic Ocean. The improvement of DAC with respect to IB over the shelves appears also in the η monthly grids, where meridionally oriented patterns of η are evidently  Table A1. Standard deviations of the three time series of along-track η , averaged over the East Siberian Sea box (Fig. A1), using uncorrected The basic cost function applied in DIVA to derive the analysis field ϕ is expressed as a sum of terms which constrain the solution as follows (Eq. B1): where d i is the observation at the location x i , µ i is the weight associated with each data point and · 2 is a norm operator. The two fundamental properties minimized by this cost function (B1) are: 1) the deviation of the analysis field from observations 660 (first term) and 2) abrupt changes in the analysis field (second term). Additional constraints J c (ϕ) can be applied.
The first term minimizes the misfits between data and analysis; weight associated to the misfits µ i are directly proportional to the signal-to-noise ratio λ. λ is to be interpreted both as a measure of the observational error and as an indication how well data represent the final analysis field (e.g., instantaneous measurements are not a good representation of a long term mean).