Geostrophic currents in the northern Nordic Seas from a combination of multi-mission satellite altimetry and ocean modeling

A deeper knowledge about geostrophic ocean surface currents in the northern Nordic Seas supports the understanding of ocean dynamics in an area affected by sea ice and rapidly changing environmental conditions. Monitoring these areas by satellite altimetry results in a fragmented and irregularly distributed data sampling and prevents the computation of homogeneous and highly resolved spatio-temporal datasets. In order to overcome this problem, an ocean model is used to fill in data when altimetry observations are missing. The present study provides a novel dataset based on a combination of along-track satellite-altimetry-derived dynamic ocean topography (DOT) elevations and simulated differential water heights (DWHs) from the Finite Element Sea ice Ocean Model (FESOM) version 1.4. This innovative dataset differs from classical assimilation methods because it substitutes altimetry data with the model output when altimetry fails or is not available. The combination approach is mainly based on a principal component analysis (PCA) after reducing both quantities by their constant and seasonal signals. In the main step, the most-dominant spatial patterns of the modeled differential water heights as provided by the PCA are linked with the temporal variability in the estimated DOT from altimetry by performing a principal component synthesis (PCS). After the combination, the annual signal obtained by altimetry and a constant offset are re-added in order to reference the final data product to the altimetry height level. Surface currents are computed by applying the geostrophic flow equations to the combined topography. The resulting final product is characterized by the spatial resolution of the ocean model around 1 km and the temporal variability in the altimetry along-track derived DOT heights. The combined DOT is compared to an independent DOT product, resulting in a positive correlation of about 80 %, to provide more detailed information about short periodic and finer spatial structures. The derived geostrophic velocity components are evaluated by in situ surface drifter observations. Summarizing all drifter observations in equally sized bins and comparing the velocity components shows good agreement in spatial patterns, magnitude and flow direction. Mean differences of 0.004 m s−1 in the zonal and 0.02 m s−1 in the meridional component are observed. A direct pointwise comparison between the combined geostrophic velocity components interpolated onto the drifter locations indicates that about 94 % of all residuals are smaller than 0.15 m s−1. The dataset is able to provide surface circulation information within the sea ice area and can be used to support a deeper comprehension of ocean currents in the northern Nordic Seas affected by rapid environmental changes in the 1995–2012 time period. The data are available at https://doi.org/10.1594/PANGAEA.900691 (Müller et al., 2019). Published by Copernicus Publications. 1766 F. L. Müller et al.: Geostrophic currents in the northern Nordic Seas

(2013)), while keeping a coarser resolution in other regions of the Earth (e.g. FESOM, Wang et al. (2014) or MPAS-Ocean model, Ringler et al. (2013)). One of the unstructured model is the Finite Element Sea ice Ocean Model (FESOM). For the northern Nordic Seas, an eddy-resolving configuration has been developed, enabling the simulation of small-scale eddies down to 1 km (Wekerle et al. (2017)). Besides total ocean current velocities including wind-driven and geostrophic components, FE-SOM includes sea surface heights with respect to the bottom ocean topography, which can be also seen as an estimation of the 5 dynamic ocean topography. Applying the gradient of these differential water elevations leads to the computation of simulated geostrophic currents. In contrast to observational based data, models show consistent spatio-temporal resolutions and enable investigations of ocean surface currents under the sea ice layer. However, they are limited to a fixed defined mathematical background and function as an assumption to the reality.
The current publication aims to present an innovative combined data product, based on the advantages of both, simulated 10 and observed datasets. In contrast to other commonly used datasets or assimilation methods, the introduced product is mainly focused on the observational side by filling in modeled DOT elevations, where altimetry data is missing or corrupted. Several investigations and consistency checks have been made by  concluding with a good agreement of simulated and observed DOT in terms of the most dominant seasonal signals and spatial patterns aiming at a combination of the temporal variability provided by altimetry along-track derived DOT elevations with simulated spatio-temporal homogeneous 15 DOT heights of the model. The combined dataset obtained is characterized by the spatial homogeneous resolution of the model and the temporal variability of altimetry derived DOT elevations. This enables further in space and time consistent studies of geostrophic surface currents in sea ice regions and may help to deepen the knowledge about polar ocean current dynamics.
The dataset is based on a combination of multi-mission satellite altimetry data from the ESA missions Envisat as well as ERS-2 and the eddy-resolving model, FESOM, covering a period of about 17 years. The combination approach is based on 20 the common known Principal Component Analysis (e.g. Jolliffe (2002), Preisendorfer (1988)), which is successfully applied in historic sea level analyses and reconstruction investigations (e.g. Ray andDouglas (2011), Church et al. (2004)).
The study area covers the Fram Strait region, the Greenland Sea and parts of the Norwegian Sea as well as the Barents Sea.
The different regions are summarized by northern Nordic Seas. In geographical coordinates the investigation area is limited to -30°W to 30°E and 72°N to 82°N. 25 The paper is structured in four sections. First, the datasets and combination method are introduced, followed by the results.
Furthermore, the combination's reliability is evaluated by comparing the obtained datasets with in-situ drifter velocities and independent satellite derived DOT products. The study closes with a summary and concluding remarks of the most significant aspects.  (right), illustrating season-dependent data coverage. a latitudinal limit of 81.5°N. The data pre-processing from ERS-2 and Envisat observed ranges to derived DOT heights and follows the descriptions of . Altimetry ranges are retracked by ALES+ (Passaro et al. (2018)) and open water/sea ice discriminated by applying the method of Müller et al. (2017). The obtained sea surface heights are reduced to DOT estimates by subtracting the high resolved and with in-situ data combined Optimal Geoid Model for Modeling Ocean Circulation (OGMOC), developed up to a harmonic degree of 2190 (Gruber and Willberg (2018)). ALES+ has been chosen as 5 an optimal retracking algorithm due to the ability for a consistent range estimation independent of the backscattering surface (open-ocean,lead and polynya). Coarse outliers are excluded from the dataset by filtering the sea surface heights on the basis of sea level anomalies (i.e. sea surface heights minus a mean sea surface) before transforming them to physical DOT heights.
A time mean inter-mission offset is removed by taking the Envisat time series as a reference within a 6-months overlap period (2003/01 -2003/06) considering only height observations from ice-free regions in the southern part of the investigation 10 area. Before introducing the altimetry DOT elevations to the further processing steps, the ellipsoid referenced observation coordinates are transformed to consider the spherical Earth representation of the model. Figure 1 shows as an example 3 days of altimetry data during the winter (March 2004) and summer (July 2006) season. In the winter, big data gaps can be noticed close to the East Greenland coast due to the presence of sea ice in contrast to summer, when most of the data is available. The second part of the combination consists of simulated differential water heights (e.g. Figure 2) with respect to the ocean bottom topography (i.e. bathymetry). The bathymetry acts as geopotential surface, which enables a linkage to the altimetry derived DOT heights (Androsov et al. (2018)). FESOM is a global multiresolution circulation model with an included sea ice component resolving the major sea ice drift patterns. The model is based on the standard set of hydrostatic primitive equations 20 in the Boussinesq approximation (Wekerle et al. (2017)) and is characterized by an unstructured triangular mesh with 47 vertical levels (Wang et al. (2014)). The horizontal resolution in the configuration used in this study reaches up to 1 km in  the Fram Strait and northern Greenland Sea area and can be described as "eddy-resolving". Furthermore, the geographical model coordinates are referenced to a spherical Earth representation with a radius of 6,367.50 km. More details of the FESOM configuration can be found in Wekerle et al. (2017). The present study uses only daily differential water heights (DWH) of the surface level covering the period 2002-2009.

Comparative datasets 5
For validation a comparison with external generated absolute dynamic topography (ADT) elevations, from ADT derived geostrophic velocity components and to geostrophic ocean velocities reduced in-situ drifter observations is performed. The ADT data including geostrophic velocity components (Pujol and Mertz (2019)), provided by the E.U. Copernicus Marine Service Information platform (CMEMS) are characterized by a daily and 1/4 degree spatial resolution and are based on multimission altimetry data. The ADT grids are created by adding temporal variable sea level anomalies to a mean dynamic topog-10 raphy and cover the complete time period of the developed datasets. However, no ADT and current data are available in sea ice areas, which limits the comparison to ice-free regions and seasons.
Further interpolated surface drifter trajectories from CMEMS (Rio and Etienne (2018)) with a 6-hour interval are used.
Following the pre-processing steps of the drifting buoys, described by Rio and Etienne (2018), all surface drifters are analyzed concerning their drogue status and local wind slippage corrections. Besides geostrophic velocities, drifter observations include 15 a-geostrophic movements (e.g. Ekman currents, Stokes drift, inertial oscillations, local wind effects, etc.). Hence, the drifter data must be corrected enabling a comparison with satellite altimetry and simulated derived geostrophic currents. Local wind corrections, also provided by CMEMS (Rio and Etienne (2018)), are directly subtracted from the drifter velocities, considering the drogue status. The Ekman current is taken from global grids providing velocities at 15m depth (drogue on) and at the surface (drogue off) level. The computation of the Ekman fields follows the explanations and processing scheme of Rio and Hernandez 20 (2003) and Rio et al. (2014). The 3-hourly available Ekman grids are downloaded from the GlobCurrent data repository and have a spatial resolution of 1/4 degrees and a global coverage. However, grid nodes north of 78.875°N are not defined, which The Stokes drift is provided globally (Rascle and Ardhuin (2013) distributed by GlobCurrent) and applied only to undrogued surface drifter data in the same way like the Ekman fields (Rio et al., 2014). Following the suggestions of Andersson et al. (2011), the Ekman and Stokes drift reduced drifter velocities are low-pass filtered by a 25-hour cutoff, two-point Butterworth 5 filter to remove tidal and inertial oscillations. Furthermore, drifters showing observations with time gaps of more than 1 day are filtered separately (Andersson et al. (2011)).
Most of the drifter buoys observations are collected in ice-free regions affected by currents (see Figure A1). Analyzing the geostrophic amplitudes and phases, the major pathway and stream velocity of the West Spitsbergen Current is clearly identified, in contrast to the East Greenland Current, which is mostly covered by sea ice. Due to high variability, most of the

Method
In order to generate a combined spatio-temporal consistent dataset based on irregular distributed altimetry observations, it is necessary to connect the along-track derived DOT estimates with a spatially consistent modeled DOT representation to fill the observation gaps. The following chapter describes briefly the combination of along-track DOT heights with the modeled water 20 level, while keeping the spatial height reference of the altimetry observations.
The combination is mainly based on a Principal Component Analysis (PCA) transferring the method of historic sea level reconstruction (e.g. Church et al. (2004), Ray and Douglas (2011)) to the present purpose. Altimetry observed along-track DOT heights represent the temporal DOT variability, whereas the spatial signal is provided by FESOM. Figure 3 highlights the interrelationship of the datasets and gives an overview over the main processing chain. The individual work steps are described 25 chronologically. The output of the processing steps are combined geostrophic currents (cGC) and dynamic ocean topography (cDOT) data representing the temporal variability of the altimetry measurements and the spatial homogeneity of the ocean model.

Data pre-processing
The input of the data production chain are along-track DOT elevations and daily simulated finite element formulated differential 30 water heights (DWH). In order to establish an equal combination basis, both datasets are treated equally. First, they are reduced by their time mean offsets and the most dominant seasonal (i.e. annual) signal ).   In a second step, the reduced FESOM grids are introduced to a PCA in order to decompose them in a linearly uncorrelated, temporal part (i.e principal components) describing the temporal evolution, and in empirical orthogonal functions (EOF) identifying most dominant spatial structures of the time series. They are sorted in a decreasing order with respect to their contribution to the total signal variance. In order to reconstruct the original signal, the principal components and the corresponding EOFs have to be multiplied and summed up. The product of one combination pair is called Mode. This inverse process of PCA is also 5 called Principal Component Synthesis (PCS). Not necessarily PCS is used to reconstruct always the full signal, however the approach can be also limited to a certain number of retaining Modes, representing a significant percentage of the total signal.
Mathematical and functional relations are explained in Jolliffe (2002). In order to determine the number of the most significant EOFs, the root mean square error (RMSE) is computed for comparing the original FESOM DWH and the reconstructed signal.
The RMSE is computed by (Barnston (1992)): where l substitutes the original FESOM DWH and r the reconstructed grids of the day t, where the overbar is computed over all grid nodes. Figure 4 shows the evolution of the temporal amount of variance and the temporal averaged RMSE with respect to the individual number of Modes. It is decided to use 50 Modes resulting in an RMSE of about 10 mm and a summed variance of more than 99%. In the following processing steps, only the spatial signals (i.e. EOFs) of FESOM are used. In contrast, the 15 principal components, describing the temporal evolution of the different modes, are neglected.

Combination
The combination step links the pre-processed along-track DOT heights with the most significant spatial pattern obtained from the PCA of the FESOM differential water heights. The processing is based on a daily temporal resolution, including 9 days of radar altimetry data for each time step. The time steps are referred to the mean of a 9-day time span (i.e. t ±4.5d ). The combined DOT heights (cDOT) can be represented by a linear combination of n combined estimated principal components 5 and the obtained EOF grids from FESOM. The functional relation of the PCS is described in Equation 2: where n corresponds to the number of significant principal components and empirical orthogonal functions. pc combi substitutes the n unknown combined principal components and EOF i (x, y) the n most dominant spatial pattern on the FESOM grid (see sec. 3.1).

10
The principal components (pc combi ) are estimated by fitting the model EOFs to the altimetry derived DOT elevations dot res .
Therefore, the EOF grids are interpolated to the observation coordinates based on NN-Interpolation resulting in along-track sampled empirical orthogonal functions (eof i (x, y)). The solution for pc combi is then given by applying the least squares method (e.g. Koch (1999)) to Equation 3: 15 where dot res (x, y, t ±4.5d ) includes all altimetry derived DOT heights within ±4.5 days and eof i (x, y) the corresponding along-track interpolated modeled EOFs. The result are n time series of combined principal components. Related detailed background information to solve for pc combi is prepared in Appendix A.
Furthermore, a Gaussian weighting, which considers uncertainties in the altimetry DOT heights due to the presence of sea ice, is introduced to the least squares process. The individual weights are defined by using external sea ice concentration from  2. Nearest-neighbor interpolation of EOF to altimetry along-track observations (coord obs ) obtaining profiled eof .

Least squares estimation (see Appendix A) of combined principal components (pc comb ) by solving Equation 3 based on
altimetry DOT observations (dot res ) and interpolated empirical orthogonal functions (eof ).

Application of Equation 2
to obtain the combined DOT (cDOT res ) dataset in the FESOM grid (coord sim ) based on pc comb (step 3) and EOF (step 1).

Data generation
In order to reconstruct the full signal and to rescale the combined heights to the altimetry height reference, the previous subtracted altimetry time mean offset and annual signal are re-added (Sect. 3.1). In the next step, combined geostrophic currents (cGC) are obtained by computing the zonal (u g ) and meridional (v g ) geostrophic velocity components at the surface, given by where g is the acceleration of gravity (9.832 m s 2 ), f = 2Ω sin φ the Coriolis force, φ the latitude and Ω the Earth's rotation rate. ∂h denotes the horizontal gradient in x and y direction of cDOT height h. The derivatives ∂h ∂y and ∂h ∂x are solved based on the finite element method (see Appendix B) which prevents further smoothing effects. Furthermore, the geostrophic absolute velocity (A g ), phase φ g and eddy kinetic energy (EKE) can be computed by applying Equation 5.
where t substitutes the velocity at a certain time and the overbar indicates the mean velocity for a defined time period (e.g. quarterly).

Datasets
The combined DOT and geostrophic current velocity fields contain DOT heights derived from satellite altimetry and simulated differential water heights from FESOM ). The dataset spans a time period from mid-May 1995 to early

Comparison with external datasets
The produced datasets are compared with independent datasets providing daily sampled DOT heights and observations of 25 surface drifter buoys. However, it must be noted that the comparison is challenging since no dataset can be used as ground truth in the whole area of study.
In order to follow a comparison with in-situ observations, the combined geostrophic components are spatio-temporally in- and from the derived combined product. Therefore, the combination procedure is applied to the drifter epochs. This is done by interpolating the estimated combined principal components linearly to the drifter times followed by a PCS (eq. 2) and a spatial nearest-neighbor interpolation to the drifter location. The result are combined DOT heights at the drifter observation time and location. In order to compare with the geostrophic drifter measurements, the cDOT heights are transformed to geostrophic velocities (sec. 3.3).    Figure 9 shows daily three-monthly averaged EKE of the combined and ADT grids within the investigation period (1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012). The EKE results are computed by subtracting three-monthly means from the daily datasets (Equation 5). The ADT appears smoother and shows big data gaps in sea ice regions in comparison with the combined results. Furthermore, the combined eddy fields show finer eddy structures within the sea ice area and close to the Greenland coast. The cGC are characterized by a higher spatial resolution and more variability in current regions.

30
The cDOT grids are evaluated against the daily and spatial averaged time series of ADT fields. Therefore, the cDOT fields are spatially interpolated to the ADT grids. Figure 10 shows the by their mean reduced temporal evolution of both datasets.
The comparison covers the full investigation period, but spatially limited to ice-free regions.

Summary and Conclusions
The current paper presents an innovative dataset based on a combination of height observations from satellite altimetry with 5 spatial information provided by an ocean model (FESOM). In case of altimetry data, an open water classification procedure is applied in order to exploit along-track water height measurements within the sea ice area. Furthermore, height offsets between the open ocean and the sea ice area are removed by using one single retracking algorithm.
The combination approach takes advantage of the principal component analysis, especially the separation of the model data into its most significant spatial patterns and temporal components with respect of the total variability. The 50 most dominant 10 patterns (EOF) are used to combine them with ERS-2 and Envisat observed along-track DOT heights in order to fill in observational gaps and to enable investigations based on a homogeneous DOT representation. In detail, the spatial information from FESOM and the temporal variability from altimetry are linked. The height level of the final product is given by altimetry by readding the previous estimated and subtracted annual signal and constant offset, since the model height reference is not clearly defined. Whereas the obtained spatial resolution is given by FESOM, characterized by local refinements in ocean current active The solution for pc combi is given by applying least squares (Koch (1999)): where β substitutes pc combi , l the observations dot res and E the interpolated empirical orthogonal functions, eof . The

Appendix B: Derivation of Finite-Elements in FESOM
The FESOM configuration that was used is based on a finite element formulation. Regarding the spatial discretization, the global ocean is discretized by using tetrahedral elements. These elements are constructed by first generating a surface triangular mesh. In the vertical, z-layers are used. The resulting vertical prisms are then cut into three tetrahedrals. In the finite element method, variables are approximated as linear combinations of a finite set of basis functions {N i }. Regarding the choice of these basis functions, FESOM uses a P1-P1 discretization, meaning that piecewise linear basis functions are employed for both sea surface height η and horizontal velocity u: where N 2D and N 3D denote the number of 2D and 3D nodes respectively. The ith basis function N i is equal to 1 at node i and linearly vanishes to 0 within elements containing this node.
Derivatives are computed by transformation to a reference element. In 2D, we consider the reference elementK defined by nodesâ 1 = (0, 0),â 2 = (1, 0) andâ 3 = (0, 1). As local 2D basis functions defined onK, we choose the first order polynomials can be mapped on the reference elementK by affine-linear transformation: with B = (a 2 − a 1 , a 3 − a 1 ) and d = a 1 . When computing the gradient of a variable φ on the reference elementK, we obtain: Thus, the gradient in the physical domain can be expressed as: We now compute the gradient of η on element K by inserting φ = 3 i=1 η i N i in the above equation: