Development of a continuous spatiotemporal finite element-based representation of the mean sea surface

The mean sea surface (MSS) has an important role, both, in the calculation of the mean dynamic topography and in the area of sea-level change as a reference surface. This paper presents a new approach to estimate a continuous spatiotemporal MSS from along-track altimetric sea surface height measurements. A parametric function continuously defined in the spatial as well as temporal domain is constructed from a C1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C^1$$\end{document}-smooth finite element space to represent the MSS. Least-squares observation equations are set up, to estimate the unknown scaling coefficients from the sea surface height measurements as collected by altimetric exact repeat missions and geodetic missions. An advantage of the proposed method is that the surface is represented by an analytic model and the unknown parameters can be physically interpreted. Whereas the static component of the function represents the MSS, the temporal component is used to absorb the ocean variability. Within this initial study, 10 years of satellite altimetry over the period 2010–2019 are analyzed in a small study region around the Agulhas Current. To obtain the best possible data coverage, all missions available via AVISO, which provide data in the study region and period, are included, i.e., geodetic missions and mission phases as well as exact repeat missions. Besides the static MSS, the temporal component, which is co-estimated to absorb the dominating ocean variability, is modeled with different basis functions to study their performance. On the one hand, global basis functions considering a linear trend and periodic functions are compared with B-Spline basis functions. The comparison of the static component to the global CNES_CLS15 MSS shows a reasonable agreement with a root-mean-square error in the range of 1–4 cm for the well-suited model configurations. To validate the modeling approach and the different analyzed configurations, the temporal model component is compared to gridded sea-level anomaly products. Although it is not (yet) a target quantity, the analysis can serve as a quality check of the MSS and the proposed modeling approach as well. It is shown that in regions with relatively low ocean variability the combination of a linear trend with an annual period is well suited to model the dominant temporal signal, whereas it is not sufficient in regions with strong ocean variability, e.g., close to the Agulhas Current. In those regions, the scenario which utilizes B-Splines in the temporal domain performs significantly better. In general, it is demonstrated that the proposed approach can be an alternative to the well-established MSS estimation procedures.


Introduction
The mean sea surface (MSS) as the spatial and temporal average of the actual sea surface is used as a reference surface in oceanography, but also in geodesy. In this context, it finds applications, for example, in the computation of bathymetric maps (see, e.g., Dixon et al. 1983;Smith and Sandwell 1994), in the geodetic computation of the mean dynamic ocean topography, which can be used to detect and analyze global long-term stable ocean currents (e.g., Knudsen et al. 2011;Becker et al. 2014;Mulet et al. 2021), or in general sea-level studies (e.g., Chen et al. 2000;Uebbing et al. 2019).
A common method of calculating the MSS in geodesy is a grid-based spatial and temporal averaging of along-track sea surface height (SSH) observations measured by satellites equipped with radar altimeters in a given period. A major challenge in calculating the MSS is the temporal variability of the sea surface observations (e.g., Andersen and Knudsen 2009;Schaeffer et al. 2012;Agha Karimi et al. 2020), which-due to temporal aliasing-can cause local biases in the MSS models.
At this point, it becomes important to distinguish between the observations of the exact repeat missions (ERMs) and the geodetic missions (GMs). Whereas the ERMs like Topex/Poseidon or the Jason family have a short repeat phase of just 10 days, the GMs have a long repeat period, which is typically above one year. Consequently, ERMs have a poor spatial coverage, as the track spacing is large. On the contrary, GMs are important for the spatial coverage, but suffer from a poor temporal resolution (e.g., Andersen et al. 2021).
For the ERMs, the elimination of the temporal ocean variability (OV) from the SSH observations for MSS modeling is straightforward. Temporal averaging of complete cycles, linked to a reference track, directly removes the periodic short-term signals and yields the mean profile of the MSS along the reference track (e.g., Andersen and Knudsen 2009;Schaeffer et al. 2012;Agha Karimi et al. 2020). The spatial processing can be performed in an independent second step, e.g., the prediction to a smoothed regular spatial grid. Unfortunately, this approach cannot be applied to the observations of GMs, as the-especially seasonal-variations do not cancel due to the repeat cycle larger than a year. Consequently, spatial and temporal averaging must be performed either i) in a joint analysis of ERM and GM data or ii) the temporal characteristics must be determined from complementary information (e.g., ERM-only or model data) and reduced within a preprocessing step ahead of the spatial analysis (e.g., Schaeffer et al. 2012;Agha Karimi et al. 2020) from the SSH measurements observed by the GMs.
Several methods have been developed to deal with the temporal variability of SSHs, which is mainly attributed to oceanic variability. In Agha Karimi et al. (2020), a spectral analysis is used to correct CryoSat-2 altimetric observations to estimate a regional MSS around Australia. Annual and semiannual amplitudes and phases are estimated from ERM data and used to reduce the ocean variability after a spatial interpolation to the CryoSat-2 tracks. The final MSS is obtained by spatially averaging the corrected data in regular 0.1 • bins.
To compute the DNSC08 (Andersen and Knudsen 2009), a two-step procedure is applied. In a first step, a longwavelength MSS is derived from the temporally averaged mean profiles of the ERMs. The short wavelengths of the MSS are derived mainly from the GM observations in a remove-compute-restore approach applied in the second step. Both are merged together using crossover analysis followed by an interpolation to a regular 1/60 • grid. For the computation of the successor MSS (DTU13 and DTU18, cf. Andersen et al. 2015), similar methods were applied. But, in Andersen et al. (2018) interpolated daily sea-level anomaly (SLA) maps (Traon et al. 1998;Ducet et al. 2000) as distributed, for instance, by AVISO are utilized to correct the GM observations for inter-annual and seasonal oceanic variability. Schaeffer et al. (2012) and Pujol et al. (2018) apply a slightly different approach to reduce the oceanic variability for determination of CNES_CLS11 and CNES_CLS15 MSS models. Schaeffer et al. (2012) uses a collocation technique, the so-called optimal analysis (OA), which enables correcting a given SSH for OV derived from other ERM altimeters that serve as reference. Finally, the corrected data are predicted to a regular grid using covariance functions for the sea surface and the measurement errors (Le Traon and Dibarboure 2004). The determination of CNES_CLS15 MSS from 20 years of satellite altimetry follows a similar approach, but instead of OA, OV is reduced by the delayed-time DUACS Level 4 gridded SLA DT2010 maps (Traon et al. 1998;Ducet et al. 2000). Jin et al. (2016) use both, corrections derived from monthly SLA maps (outside ±66 • latitude) and corrections derived from ERMs (inside ±66 • latitude), to correct the observations collected by the GMs. After a crossover adjustment, the gridded MSS is computed by least-squares collocation.
The goal of this manuscript is to present a proof-ofconcept study to estimate a continuous model of the MSS from altimetric SSH observations using time-variable C 1smooth (i.d. one times continuously differentiable) finite elements as basis functions. Instead of eliminating the temporal OV in a preprocessing step, the temporal signal is parameterized and co-estimated. SSH observations taken by ERMs and GMs are jointly analyzed within a single adjustment. Whereas static parameters describe the MSS at a specific reference epoch, the time-dependent parameters represent a filtered model of the temporal OV. Pre-defined sets of temporal basis functions are combined with the spatial finite elements, to construct various spatiotemporal models with different filter characteristics in the spatial as well as the temporal domain. As the analysis can be represented in a one-step least-squares estimation, a stochastic model of the SSH measurements can be rigorously propagated to the MSS parameters in a straightforward manner.
The manuscript is organized as follows: In Sect. 2, the theoretical background to construct the C 1 -smooth finite element space is provided. Based on the finite element space, the least-squares observation equations of a static approximation procedure are constructed. These are extended by a deterministic model which captures the temporal variations and introduces the new parameters of the spatiotemporal model. Section 3 introduces the study region and the data sets analyzed within the different estimation scenarios. Different solutions for the MSS and estimates for the OV are derived and discussed. In Sect. 4, some additional analyses are presented to highlight the potential but also the limits of the proposed procedure and the setup chosen in this initial study. This paper ends with a summary and provides some conclusions in Sect. 5. Additionally, some aspects for future studies and refinements are presented.

Theoretical background
As discussed in Sect. 1, the goal in this study is to represent the MSS by a continuous mathematical function. There is no accessible closed expression that is directly derived from physical laws. Instead some function space defined over the domain of interest must be constructed, from which a concrete function can then be estimated in a least-squares adjustment.
To include temporal variability of the sea surface in the model, two independent function spaces, one defined in the spatial domain and the other in the time domain, are combined via the tensor product of their respective bases.

Finite element spaces for spatial approximation tasks
In consideration of the requirements that come with the ocean surface estimation, namely coastal boundaries and relatively high spatial variability, a sensible choice is finite elements on a triangulation in R 2 . For constructing these so-called finite element spaces (FESs), the entire domain of interest is partitioned into a finite number of (in this case triangular) subregions, each of which has its own locally defined basis functions and corresponding degrees of freedom (also called parameters). The complete definition of a finite element is then given by the shape, the local degrees of freedom and the local function space spanned by the basis functions of such a subregion (following Ciarlet 1978). By means of an appropriate choice of local basis functions, these are usually taken from the literature, a piecewise function space with at least C 0 -continuity is achieved. Higher orders of continuity are possible and depend on the finite elements actually used. Given a FES definition, a function can then simply be written as where i ∈ I s describes the indexing of all I s piecewise defined basis functions b s,i : R 2 → R, and a s,i ∈ R the corresponding degrees of freedom/scaling coefficients for modeling a spatial signal s. The set {b s,i } of all basis functions is the basis of the FES. A relational measure for the spatial resolution within a region of fixed size is the number of degrees of freedom, equivalent to the dimension of the finite element (sub-)space over that region. Adopting the spatial variability of the expected signal is in general done on a global scale (meaning over the entire domain) by selection of some specific finite element and then refined on local scales by adjusting the size of subregions.

Selected finite elements applied in this study
In this study, we consider local basis functions on triangles, or finite elements, which yield a C 1 -continuous FES, i.e., the first derivative of the function in the entire domain is continuous. The choice of higher continuity comes from the reasonable assumption that the ocean surface is smooth on the scales observable with satellite techniques, in this study specifically satellite altimetry.
Although different elements are implemented, focus in this study is on the Argyris element (cf. Argyris et al. 1968). It is the element with lowest degrees of freedom which guarantees C 1 -continuity while spanning a complete polynomial space. In particular, this is the local space of polynomials of degree 5 with 21 degrees of freedom including the function value, the two first and the three second derivatives for each of the three nodes, as well as the three normal derivatives on the centers of the edges. Thus, in addition to the general C 1continuity, this element guarantees continuity of the second derivatives in the nodes of the triangulation.
In comparison, the class of macro-elements exist, where the previously described technique to construct a full FES is similarly applied (one might call it recursively) to an individual element. One example is the Hsieh-Clough-Tocher (HCT) element (Clough and Tocher 1966), which provides the same continuity with only 12 degrees of freedom. Based on a further sub-partitioning of a triangle into three subtriangles, each of which is domain to a polynomial of degree 3, a C 1 -continuous piecewise local function space is defined by implicit constraints in the 12 basis functions. These correspond to the 12 degrees of freedom as defined for the whole element, namely the function value and the two first derivatives in each node, as well as the normal derivative at the center point of each edge.

Triangulation
The second major defining factor of a FES is the partitioning of the domain of interest into a mesh, here by means of triangulation. It has an immediate effect on the achieved spatial filtering via the size of individual triangles. The coarser the mesh in a particular region, the more the signal is filtered and vice versa.
A sophisticated method to generate meshes is developed and implemented in the software library jigsaw (Engwirda 2014(Engwirda , 2017(Engwirda , 2019 and can, for example, be used through its bindings to scripting languages (matlab, octave or python). The software uses a Frontal-Delaunay refinement technique which allows for the generation of high-quality unstructured spheroidal Delaunay triangulations and uses a hill-climbing-type optimization technique for further grid quality improvement. It brings the desirable properties of being automatic in its generation of meshes (given some configuration input) and enabling easy local adjustments of triangle density/size as well as geometric constraints. In this proof-of-concept study, planar meshes in plate carrée projection with constant mesh spacing for the whole domain are used. The study regions are restricted by polygons derived from coastline data and straight boundaries over the ocean.

One-dimensional models for the time domain
For modeling the temporal behavior of the sea surface, we make the assumption that changes in time are continuous and separable from spatial variability. The latter means that at any fixed point in space, the OV can be represented by a function with the only variable being time (similar to Jin et al. 2016; Agha Karimi et al. 2020). As such, we consider a separate function over time domain where d s, j : R → R are the J s basis functions and c s, j ∈ R once again the corresponding degrees of freedom/coefficients. We consequently say that {d s, j } is the set of all J s temporal basis functions.

Temporal basis functions with global support
It is necessary to identify basis functions that can model the dominant expected and resolvable temporal variations. A first set of candidate basis functions follows from the studies by Jin et al. (2016) and Agha Karimi et al. (2020) and is a temporal model which covers a linear trend and harmonic terms for the annual (and semiannual) signal. As the assumption of a strictly linear trend is only sensible for sufficiently short periods of time, a simple transformation is introduced that maps a point in time t with respect to some reference epoch t 0 , i.e., Δt = t − t 0 . This is analogous to a linearization at t 0 .
The set of temporal basis functions is then given by d s, j := {Δt, cos(ω 1 Δt), sin(ω 1 Δt), . . . , cos(ω L Δt), where ω l = 2π T l (l ∈ L) is the angular frequency and L describes the set of harmonics with the period length T l to be estimated. Since all selected basis functions have zero mean (assuming a sampling of full periods and a choice of t 0 in the center of the time series), it is guaranteed that this kind of basis functions does not absorb static signal.

Temporal basis functions with finite support
In contrast to the global basis functions (cf. Sect. 2.2.1) which are rigid, basis functions with local support, such as the finite elements in spatial domain, allow a more flexible approximation. Thus, other feasible choices for the temporal basis functions are one-dimensional finite elements, i.e., piecewise polynomials, or spline function spaces over R, like B-Splines (see for instance De Boor 2001;Fahrmeir et al. 2021).
Equivalent to (3), the set of basis functions for B-Splines can be written as where B p k (Δt) describes the k-th B-Spline (k ∈ κ) of degree p of the nodes in the temporal domain κ (e.g., Fahrmeir et al. 2021, Chap. 8). In this study, uniform B-Splines of degree p = 3 with equal node spacing of Δκ in the entire temporal domain are used. Next to the higher flexibility another advantage of B-Splines is the well-known low-pass filter characteristic which only depend on the node spacing and their order (cf. Sünkel 1985;Köhler et al. 2019). Advantages are offset by significantly more degrees of freedom and a reduced physically interpretability of the coefficients. Furthermore, as the mean of the basis functions is not zero, additional constraints are required (cf. Sect. 2.6) to avoid the leakage of static signal into the spatiotemporal model.

A separately continuous spatiotemporal model
The spatiotemporal model space is then composed by the tensor product of the spatial and temporal model spaces. Given the two sets of basis functions {b s,i } and {d s, j }, with the tensor product for real-valued functions defined as the combined set of the I s J s basis functions to model a signal s is In terms of the concrete basis functions {d s, j }, it is for the case of global temporal basis functions (cf. with the degrees of freedom/coefficients e s,i, j . In a more explicit form with function arguments and using (7) (indicated by the superscript g), this becomes or again equivalent for B-Splines (indicated by superscript l) using (9)

A model function for MSS estimation from SSHs
To model the MSS by a continuous function, a FES is constructed and tailored to the expected MSS signals and the distribution of observations. As there is no time dependency, (1) is used to obtain the model function for the MSS, i.e., To absorb the temporal OV, a spatiotemporal model function using (9), which either follows from (10) or (11), is used and represented as the generic function Here, the FES is tailored to the observational geometry of the ERMs to compensate the OV and to avoid leakage into the MSS. In case of a zero mean function f OV , the SSH can be then approximated by the sum Here, f SSH contains the spatiotemporal model component to reduce the ocean variability, especially from the SSH observations collected by the GMs.

Estimation approach
To estimate a MSS model from stochastic SSH observations L, here in terms of instantaneous satellite-based altimetry, the standard least-squares adjustment method is applied (e.g., Koch 1999, Chap. 3).
In this initial study, the SSH observations are assumed to be independent and identically distributed with a variance of σ 2 0 , i.e., the covariance matrix is The realized measurements l are further assumed to be point-wise. The p-th measurement is taken at location with N the number of altimetric SSH observations. As will be further discussed in the following section, data from multiple satellite missions are used to estimate the models. Although a priori intermission bias corrections are applied by the data provider (cf. Sect. 3.1), the observation equations include additional bias correction parameters o m for the different missions, to account for a local refinement. A single mission with a fixed bias of 0 m is selected as reference.
The scenarios S * (cf. Sect. 3.2) do not account for the temporal model component. The least-squares observation equations for an observation l p taken by mission m simply read (cf. (1)) where v p is the least-squares residual. For the scenarios A * , * , A * , * and B * , * , the least-squares observation equations derive with equation (14) as The vector of unknowns x is then composed of the coefficients, i.e., {a MSS,k , e OV,i, j , o m } for (17) or {a MSS,k , o m } for (16), respectively, for the spatiotemporal or spatial models and is estimated in a linear Gauss-Markov model.

Constraints for the B-spline model
The temporal modeling with the B-Splines increases the flexibility and enables the modeling of small-scale OV. But to separate the SSH observations into the MSS signal and the temporal ocean variability, it has to be considered that the B-Spline function, in contrast to the global temporal basis functions cf. (10), can also absorb static signal.
To prevent this leakage, additional constraints to the spatiotemporal model coefficients are required to obtain a zero mean. Thus, the mean of the function in the temporal domain 1 T T f l OV (λ, φ, Δt)dt has to be zero. From this postulation, the linear restrictions to the parameters k∈κ e OV,i,k+1 = 0 ∀i ∈ I OV (18) can be derived, see Appendix A for details.
To increase the numerical stability, it is beneficial to constrain the temporal start node κ 0 and end node κ m of the B-Spline functions. The second derivatives of the first and last nodes can be set to zeros (cf. Köhler et al. 2019), so that the 2 × I OV linear constraints follow for all start and end nodes; see Appendix B for details. Both sets of linear restrictions are applied in the parameter estimation of the scaling coefficients in a Gauss-Markov model with constraints (e.g., Koch 1999, Sect. 3.2.7).

Regularization strategies
Although basis smoothness characteristics of the spatial functions g s (x, y) follow from the definition of the FES, additional stabilization and smoothness can be achieved applying regularization strategies such as Tikhonov-like regularization (cf. Tikhonov et al. 1977) which can be used to regularize (local) parameters directly. Regularization at the level of linear functionals of g s is used in the context of approximation, such that in addition to minimizing the weighted squared sum of residuals, the weighted norm of the spatiotemporal scaling coefficients i∈I s j∈J s e 2 s,i, j is minimized. In the context of MSS approximation using FES, applying regularization strategies helps to solve numerical issues due to individual poorly determined parameters (especially in coastal and boundary regions) or spatial and temporal data gaps.

Data and study area
To obtain a well-suited model of the MSS, it is important to have data with, both, high temporal and high spatial resolution (cf. Sect. 1). Thus, altimetric SSH measurements from all satellite missions available at AVISO are selected. The time period considered in this study is 2010 to 2019 inclusive. The high temporal resolution (from ERM data) is realized using the measurements of the regular mission phases of Jason-2, Jason-3, as well as the Jason-1 and Jason-2 interleaved mission data, Envisat Nominal and Extension Phase, SARAL/AltiKa, HY-2A and Sentinel-3 A and B (see, e.g., ESA 2020b, c, d). Besides the Jason-1, Jason-2 and HY-2A GM phases which fall into the period of interest, (pseudo) Low Resolution Mode CryoSat-2 and SARAL Drifting Phase observations are used for their high spatial resolution (see, e.g., ESA 2020a; Andersen et al. 2021). Due to the nearly complete temporal coverage of Jason-2 in the analysis period and its high temporal resolution, it is used as the reference mission for bias estimation. Observations from the L2P data product as processed on behalf of CNES (Centre National d'Etudes Spatiales) SALP project and distributed by AVISO+ 1 were used. The applied instrumental and geophysical corrections such as for orbit, intermission bias or tidal errors can be found in the AVISO L2P SLA product handbook (see AVISO 2017).
Some of the mission characteristics are shown in Table 1 (cf. Andersen et al. 2021;ESA 2020b, c, d, a). Whereas all missions have similar along-track resolution of about 7 km as 1 Hz SSH data are used, the ground track spacing significantly differs from 8 to 315 km (cf. Table 1).
To keep the huge amount of along-track observations in the selected time period manageable, a small study area is selected to test and analyze the performance of the proposed approach. This area (see Fig. 1, blue polygon), the South Atlantic Ocean and Indian Ocean south of Africa with the Agulhas Current, is selected because of its individual characteristics.
Although having a small extend, it contains regions with low as well as high temporal ocean variability, especially due to the Agulhas Current. This allows for testing the approach in various configurations, e.g., the estimation of models on finer meshes, with limited computing resources. The number of used observations inside the study area is shown in the last column of Table 1 for each mission. In total, there are more than 9 × 10 6 observations. These observations are  Fig. 1 The study area south of Africa (blue polygon) and an exemplary mesh with a target edge length of 175 km with the CNES_CLS15 MSS (a) and the individual subregions R 1 -R 4 and points P 1 -P 4 (b) for more specific analyses almost equally distributed over the time period and within a year. Thus, neither a year nor a month is overrepresented in the data. Figure 1 shows an exemplary mesh with a constant target edge length of about 175 km generated for the study area (Fig. 1a) and in blue the boundary used to compute the triangulations (Fig. 1b,

Results for different scenarios
With the proposed approach presented above, it becomes possible to estimate either a static MSS model (Eq. (16)) or a MSS model which co-estimates the temporal model component (Eq. (17)). This study puts the focus on the parameters which represent the static MSS using the specific models of Eq. (10) and (11). The reference epoch t 0 is chosen as January 1, 2015, which is in the middle of the study period.
The numerical results obtained in this study are evaluated (cf. Eq. (12)) on the grid as provided by the CNES_CLS15 MSS (Pujol et al. 2018) and compared to that static MSS product converted to the same reference epoch.  (13)).
To demonstrate the flexibility of the proposed approach, seven different estimation scenarios with different spatial and temporal characteristics are studied (cf. Table 2). Whereas the S models are estimated without the temporal model extension (Eq. (16)) as a kind of baseline, the A models are estimated with a linear trend and one or more given periods (T 1 describes the annual period and indicates T 2 annual and sub-annual periods, Eq. (17) and (10)). The last model B 25,175 is estimated using B-Splines for the temporal model (Eq. (17) and (11)) with a constant node spacing of Δκ. It serves as an outlook to demonstrate the flexibility and the potential of the proposed modeling approach. In this study, we only present results which use the Argyris element for spatial modeling. The alternative HCT element, mentioned in Sect. 2.1, is an interesting alternative for future studies.

Construction of the finite element spaces
Meshes for two different FES which represent different spatial resolutions are constructed with jigsaw and analyzed in the following sections with respect to its applicability to represent a) the static MSS signal and b) the ocean variability. The first fine mesh has a target edge length of 25 km and the second coarse mesh about 175 km in the study region (cf. Fig. 1a).

Baseline scenarios without co-estimating the ocean variability
The initial scenario S 25 (Fig. 2a) is estimated on the fine mesh without a temporal model component (Eq. (16)). The result is used to demonstrate the general applicability of the FES to model the MSS signal and that the observation sampling is sufficient to derive a stable solution. Furthermore, it serves as a baseline scenario for the models which coestimate the temporal ocean variability. Figure 2a shows the difference between S 25 and the reference MSS. The differences are in the range of ±10 cm, and no obvious MSS residual signal is visible. Lowest differences occur around the repeat tracks of the ERMs.
Different statistics of the differences can be seen in Table  3. Next to the median, the RMS of the entire study region as well as the subregions R 1 to R 4 , as defined in Fig. 1b, is provided. Scenario S 25 shows an overall RMS of 3.3 cm and a better agreement in regions R 1 and R 2 , which are of lower temporal variability compared to the regions R 3 and R 4 of higher variability.
This shows that the Argyris element on the 25 km mesh is in general well suited to model the static MSS (i.e., centimeter level as visible in region R 1 ) in general, but a better modeling of the ocean variability is required. To illustrate the impact of using the coarser mesh which is tailored to the cross-track resolution of the ERMs, scenario S 175 is estimated on the 175 km mesh and still without a temporal model component. Table 3 shows that the RMS of the differences of nearly all regions increased by a factor of three, which indicates that the spatial resolution is insufficient to model the MSS signal.

Simple FES co-estimating the ocean variability
The most straightforward choice for the FES for the temporal model component (cf. Eq. (13)) is the same as used for the static MSS. Consequently, Fig. 2b and 2c shows the MSS results for scenarios A 25,25 and A 175,175 where linear trends and annual periods T 1 are co-estimated on the same mesh to absorb the temporal variability.
The differences of both scenarios S 25 (Fig. 2a) and A 25,25 (Fig. 2b) show the same general structures, but the difference for A 25,25 seems more noisy, which is supported by a slightly larger RMS (cf. Table 3). The static component requires the high spatial resolution of the 25 km mesh to model the high-resolution, mostly geoid, signal. But, the spatial resolution seems too high for a stable estimation of the temporal component, which is limited by spatial sampling of the ERM. The observations by the ERM cannot support the temporal finite elements inside the diamond structures. This causes aliasing because of the poor temporal coverage of the GM observations which results in a poor overall performance. Figure 2c shows the differences of A 175,175 to the reference model. Although the differences have a smoother character, the order of magnitude becomes larger, as for the S 175 scenario, the 175 km mesh is insufficient for the static signal. The RMS of this scenario is slightly worse than the static scenario on the same mesh. But, as will be demonstrated in Sect. 3.2.3, the temporal model shows a better agreement with reference data.
These results indicate that due to the wider ground track spacing of the ERM, it is feasible to combine different spatial resolutions (i.e., meshes with different triangle sizes). A high spatial resolution is required for the MSS and a lower spatial resolution for the temporal signal to prevent aliasing due to the temporal sampling of the GM and to avoid over-parameterization. Thus, the next sections combine the advantages of both FES by using the fine and coarse FES for the MSS and the temporal variability, respectively.

Combination of finite element spaces
The scenarios A 25,175 andÂ 25,175 realize the higher spatial resolution to model the MSS using the fine 25 km mesh in Eq. (12) and the lower resolution to model the ocean variability using the coarse 175 km mesh in Eq. (13). Figure 2d shows the differences of the MSS derived in scenario A 25,175 . This composition of the different spatial and temporal resolution reduces the RMS to nearly a third with respect to A 175,175 and about 7mm compared to A 25,25 for the entire study area. Compared to the static S 25 scenario, the co-estimation of the temporal model improves the consistency to the reference MSS by 4mm which is more than 10%. As supported in Table 3, the reduction is mainly visible in the regions R 3 and R 4 with strong ocean variability where the improvement is even larger. This shows that modeling the trend and annual period helps to reduce the temporal aliasing and consistently improves the RMS of the differences to the reference model. But, large differences remain in areas with strong variability, indicating that the temporal modeling is not yet sufficient. Therefore, the scenarioÂ 25,175 extends the modeling of the annual period by additional sub-annual periods (cf. Table 2) to increase the flexibility. The comparison of Fig. 2e and 2d shows again visual reduction of the differences due to the additional periods. This is confirmed by the RMS over the entire study area, which improves again by 2 mm (1-2 mm in all subregions). But, the visible rough differences especially in the area of the Agulhas Current show that the temporal model might not yet be sufficient to completely reduce the influence of the temporal variability.
With this kind of modeling, the estimated MSS reaches a consistency with the reference model of 1-2 cm in regions of low ocean variability and 3-4 cm in regions of high ocean variability. This is good considering (a) different data sets from different periods were used and the adjustment of the reference epoch is not perfect and (b) the CNES_CLS15 MSS uses 20 years of observations compared to the ten used in this study. To quantify an expected consistency, the RMS between the CNES_CLS15 MSS model and the alternative DTU15 MSS is 2.8 cm in the study region, which supports the obtained consistency level of 2.7-2.9 cm of the A scenarios (similar to the comparison results shown in Pujol et al. 2018).

Analysis of the temporal parameters
A unique characteristic of the presented method is that due to the chosen parameterization of the finite elements and the temporal model, the parameters for the trends and the periodic functions are interpretable. As the parameters in the static case correspond to function values, first derivatives, etc., in the vertices of the mesh, this holds true for the temporal parameters as well.
Use of Equation (13), more specifically the spatially interpolated parameters, which are e OV,i, j b OV,i (λ, φ), enable one to derive and display temporal parameters like trend and annual amplitudes and phases in the entire domain. These are shown in Fig. 3 for the scenarios A 25,25 , A 175,175 and A 25,175 , where each of the aforementioned parameters was evaluated on a regular grid. Particularly the parameters of the fine model (first row of Fig. 3) show high spatial noise. This confirms the assumption of an over-parameterization of the temporal signal in scenario A 25,25 . Consequently, data inconsistencies, static signal and noise are absorbed by the temporal model and degrade the quality of the MSS. The parameters shown in the plots of the second row are derived from scenario A 175,175 . They show smooth spatial structures as expected, which confirms that the coarse mesh is better suited for the temporal model component. The third row of Fig. 3 shows the parameters from scenario A 25,175 where static and time variable parts are estimated on different meshes. These parameters show even smoother spatial structures as for the A 175,175 , which demonstrates that the combination of the FES improves both model components, i.e., the static MSS as well as the temporal component.
One way to evaluate these results is to estimate the trend and the amplitude and phase of an annual period from the DUACS Level 4 gridded SLA DT2018 maps individually for every grid cell in a least-squares approximation. Plots of these are shown as reference in the bottom row of Fig. 3 and indicate a good consistency to the results derived from A 25,175 . A visual comparison shows a good agreement between the Fig. 3 Spatial representation of the physical properties trend, amplitude and phase displacement of the three different models A 25,25 , A 175,175 and A 25,175 (first to third row) on the grid defined by the DUACS SLA maps. As comparison, trend, amplitude and phase maps as derived from the daily DUACS SLA maps are shown in the bottom row 123 parameters and the reference plots, especially for the coarse temporal models. The same structures are visible, and the orders of magnitude are the same, although the spatial structures of the reference plots are still smoother.
The comparison of the parameters with the reference derived from gridded SLA products shows that the trend and annual signal obtained from the observations are modeled successfully in the proposed combined estimation approach. Without over-interpretation of the results, this additionally validates the approach and the derived MSS and confirms its potential. Static and long-term temporal signals have been successfully separated.

Increasing the temporal flexibility with B-splines
The remaining differences (cf. Fig. 2d (11)) which are independent of pre-defined frequencies.
In this study, uniform B-Splines of order 3 are used with a constant temporal node spacing Δκ and combined to a spatiotemporal model (Eq. (14)). Within the parameter estimation, the restrictions discussed in Sect. 2.6 are applied. To improve the stability due to the high number of parameters and the finite base of the temporal model, the regularization as discussed in Sect. 2.7 is applied. To keep the influence of the regularization small, the variance factor for each parameter group i is set to σ 2 i = 1 × 10 3 . The node spacing Δκ = 30.44d is fixed, mainly for computational restrictions and corresponds roughly to the sub-cycle of CryoSat-2. Still, the number of unknown parameters in the B-Spline scenario B 25,175 is large, i.e., 293 843, triple to quadruple that of A 25,175 andÂ 25,175 (see Table 2). Figure 2f shows the differences to the CNES_CLS15 MSS. It can be seen that the highest differences are still in the region of the Agulhas Current but comparatively smaller in contrast to the differences shown in Fig. 2d and 2e. This visual interpretation is supported by a smaller RMS of 2.0 cm for the entire study area, representing an additional improvement of 25%.
The same general features remain visible as for all difference plots in Fig. 2 and thus 2f: a larger consistency in the northern area, where less ocean variability is expected (cf. Fig. 3), and larger differences in the southern area within the Agulhas Current. As the S 175 and A 175,175 have an insuffi-cient spatial resolution for the MSS, and A 25,25 suffers from the over-parameterization, these are excluded from the discussion. Since results for A 25,175 andÂ 25,175 are very similar, conclusions are valid for both of them.
As mentioned before, the greatest consistency to the CNES_CLS15 mss is visible in region R 1 . The lowest RMS of 8mm for the entire study area is observed for scenario B 25,175 compared to 10mm for A 25,175 and 13mm for S 25 , which does not co-estimate the variability. In total, the coestimation yields an RMS reduction of 38%, which indicates that a carful modeling of the ocean variability leads to better results even in regions of low variability. Although at a higher level, the relative improvements for R 2 , R 3 and R 4 are the same: about 15% for A 25,175 and above 35% for B 25,175 compared to S 25 .

Flexibility of the model: temporal model component as time series
Within the previous sections, the primary view on the estimated model was on the static MSS in the spatial domain.
To demonstrate the flexibility of the continuous spatiotemporal model and to further validate the estimated models, the spatiotemporal models of the ocean variability (cf. Equation (13)) with its estimated parameters can be seen as a one-dimensional time series for a constant location, i.e., This time series represents the estimated temporal variability, which can be compared to time series extracted from the mean reduced DUACS SLA maps. The product provides ocean variability with a high temporal resolution of 1day. For this purpose, the temporal component of estimation scenarios A 25,175 ,Â 25,175 and B 25,175 is evaluated at different locations (cf. Fig. 1b) of the provided SLA grid and compared to the corresponding time series extracted from the DUACS product in Fig. 4.
The reference SLA (yellow) and the estimated ocean variability time series of the scenarios A 25,175 (blue),Â 25,175 (red) and B 25,175 (green) are shown for four different locations P 1 , P 2 , P 3 and P 4 with different properties, and all of them are close to the corresponding region R i (cf. Fig. 2). The locations shown in Fig. 4a and 4b are both located in and near the smoother regions R 1 and R 2 , respectively, with lower temporal variability. The SLA in Fig. 4a shows that the time series is dominated by the annual period especially in the first years, which is well represented by the estimated models of ocean variability. The time series of scenarioÂ 25,175 emphasizes that the annual period is dominant at this location. Nevertheless, the flexible scenario B 25,175 shows the best fit. Because of the small magnitude of the SLA and the dominant annual period, the RMS of A 25,175 andÂ 25,175 only slightly differs from the RMS of B 25,175 which is 1.4 cm.  Figure 4b and especially 4d demonstrates the strength of the B-Spline approach and highlights the low-pass filter characteristic of the utilized B-Splines. It is shown that nonperiodic signals can be approximated quite well, and the RMS for B 25,175 is 8.6 cm for P 4 and 7.0 cm for P 2 . Nonetheless, due to the filtering it is not yet a perfect model to represent the SLA. Consequently, we prefer to call the temporal model component a compensation model for (parts) of the ocean variability and not yet a model for SLAs. On the contrary, it becomes visible for both, A 25,175 as well asÂ 25,175 that trend and selected harmonics are not sufficient to model the small-scale variations. This is emphasized by the high RMS of 21.1 cm for A 25,175 and 20.7 cm forÂ 25,175 for P 4 and 15.1 cm for P 2 for both models. The additional periods used inÂ 25,175 only show a tiny improvement, which emphasizes that the B-Splines are preferable. P 3 is located in R 3 and is characterized by the highest variability in magnitude and includes many non-periodic signals. Supported by a spectral analysis of the SLA the annual period is not the dominating one. Due to the high variability, the amplitudes of A 25,175 andÂ 25,175 seem to be biased toward some strong high-frequency signals. High-frequency variations leak into the static model (see again R 3 and R 4 in Fig. 2d and 2e) and the linear trend. Consequently, the modeled trend and seasonal signal applied do not fit the SLA time series very well and only a minor part of the variations are captured. Even the B-Spline model B 25,175 cannot perfectly capture the high-frequency signals (controlled by the 123 node spacing). The non-perfect node spacing is also visible in Fig. 4d where the temporal resolution is not high enough to approximate the small-scale peaks.
For this first proof of concept study, our results show an impressive agreement of B 25,175 and the gridded SLA product, because the model is able to compensate the dominating signal of the ocean variability. This serves as the motivation for further studies with B-Spline models to capture the entire SLA signal by a continuous spatiotemporal model as an additional output.

Summary, conclusions and outlook
In this paper, an approach for the estimation of a continuous spatiotemporal model of the mean sea surface from altimetric sea surface height observations is proposed. The continuous model is based on C 1 -continuous finite element basis functions, which represent the spatial signal. Finite elements defined on triangulations are used (Argyris element) to assemble a continuous surface in the study area of interest. Substitution of the (unknown) scaling coefficients by one-dimensional parametric functions in the time domain constructs the proposed continuous spatiotemporal model. Within that model, the unknown scaling coefficients result from the Cartesian product of the temporal and spatial scaling coefficients of the basis functions. These are determined from the SSH observations within a least-squares parameter estimation process (cf. Sect. 2).
To demonstrate applicability, and to study the performance, of the proposed approach, it is used to estimate a continuous model of the MSS while co-estimating the spatiotemporal ocean variability. For this purpose, 10 years of satellite altimetry data collected by GMs (CryoSat-2, Jason-1, Jason-2, SARAL, HY-2A) and ERMs (Jason-1 to Jason-3, Envisat, SARAL, HY-2A, Sentinel-3A and Sentinel-3B) in the period 2010 to 2019 are used. It is demonstrated that different FESs are required to successfully separate the static MSS signal from the spatiotemporal ocean variability. Whereas a high-resolution mesh (about 25 km) is required to resolve the MSS signal, a coarser resolution mesh (about 175 km) tailored to the spatial sampling of the ERMs is required to capture the ocean variability. This avoids spatial over-parameterization of the ocean variability (cf. Sect. 3).
The proposed approach estimates a static or mean component as a C 1 -continuous finite element model and introduces the possibility to co-estimate the ocean variability as an additional spatiotemporal model component. For this component, separable functions for the temporal and spatial domain are chosen, allowing one to study different functions for the temporal domain. Initially, motivated by other common approaches, linear trends and (sub-)annual periods are considered in the form of global basis functions. Although this reduces the aliasing of temporal signal to the MSS, it is shown that a refinement of the temporal model using more flexible B-Splines as finite basis functions improves the MSS solution significantly, especially in regions of high ocean dynamics.
The temporal model component mostly absorbs long-term temporal variability, especially from the SSH measurements collected by GMs. The primary output is the MSS, which is compared to the established CNES_CLS15 MSS. After adjusting the reference epoch, a consistency of 2.7-2.9 cm RMS is achieved for the A-models. This improves to 2.0 cm for the B-Spline scenario B 25,175 . This is a good agreement considering the consistency of other MSS models (which is about 2.8 cm RMS between CNES_CLS15 MSS und DTU18 MSS in the study area) and the differences in the used data sets.
Although a good consistency to the established procedures is found, it is not possible to conclude which of the MSS models perform better or are more accurate. As the utilized altimetry data sets, selected time periods and study areas differ, no fair conclusions can be drawn. Thus, we cannot conclude that the new approach outperforms the established methods. But, from a theoretical and methodical point of view there are some advantages of the one-step procedure introduced here: -No additional preprocessing, i.e., reduction of the ocean variability, is required for MSS estimation. It is derived from the data itself in the adjustment procedure such that quality measures like covariances can be derived.
Established approaches have to (implicitly) rely on a homogeneous quality assumptions of the ocean variability products. -ERM and GM data are jointly analyzed, accounting for the acquisition time. Consequently, the reference epoch is well defined, even in case of irregular data availability and quality. On the contrary, the grid-based approaches tend to introduce biases in those cases. -Due to the parametric nature of the model, it can be synthesized on arbitrary grids and, thus, directly comes along with a tailored interpolator (in space and time) including the possibility of variance propagation. -A re-parameterization of the finite elements is on the one hand used to realize the C 1 -smoothness for the composed function. Furthermore, this results in physically interpretable parameters, including the MSS and its derivatives for the static model component, and sealevel trends, amplitudes and phases (and its derivatives) for the temporal model; in the case global temporal basis functions are used. -As the entire problem is formulated as a one-step least-squares adjustment, stochastic characteristics of the observations can be easily included and propagated to uncertainties of the estimated models. Variance propaga-tion to either the MSS, the temporal model component (OV) or the combined SSH is easily possible. Particularly for regional geodetic mean dynamic topography estimation, a fully populated covariance matrix of the MSS would be beneficial. Such model uncertainties are omitted in this publication because of the low informative value due to the insufficient stochastic model.
This work demonstrated that the approach is well suited as an alternative for MSS estimation and has the advantage of co-estimating the ocean variability guaranteeing a best possible consistency. Thus, it is less dependent on prior information and the parametric nature seems well suited to combine data sets with different quality and sampling characteristics. Our results indicate that an extension toward a continuous spatiotemporal model of sea-level anomalies becomes possible. Necessary future refinements to derive competitive and even more reliable estimates of the MSS and SLA models include: -To achieve an optimal performance, B-Splines are required for the temporal model. Due to their local support, this significantly increases the parameter space. Computational and numerical challenges arise when increasing the region and/or period of interest, which have to be resolved in future research with tailored implementations on massive parallel compute clusters exploiting the sparsity of the finite elements. -As indicated by the results, the temporal resolution defined by the node spacing is not yet optimal and, thus, not sufficient to describe small temporal variations. Nevertheless, higher temporal resolutions are possible. Theoretically, this is limited by the temporal sampling of the ERM of about 10 days. When reducing the B-Spline node spacing, i.e., increasing the temporal resolution, numerical issues due to temporal and spatial data gaps can occur. These can be resolved by smoothness constraints or more advanced regularization strategies, which have to be carefully analyzed in real data scenarios to avoid strongly biased estimates. -Another possible improvement arises from data adaptive meshes to reduce the parameter space, the numerical complexity and the risk of over-parameterization. Thus, in regions with a smooth spatial signal, the mesh can be coarser than in regions with a high spatial variability. -To obtain a meaningful uncertainty description of the estimated model, a data-adaptive stochastic model for the SSH observations is required. The least-squares residuals can be iteratively analyzed to estimated the stochastic information along the orbits for the individual altimetry missions and mission phases. This information can subsequently be used to model covariance matrices of the observations. This will again increase the numeri-cal complexity but also leads to a more realistic quality description of the MSS or SLA. -Due to its generic nature and implementation, the proposed method can be applied to different spatiotemporal surface approximation tasks.