**Evaluation and Application of Newly Designed Finite Volume Coastal Model FESOM-C, Effect of Variable Resolution in the Southeastern North Sea**

**Ivan Kuznetsov 1,\*, Alexey Androsov 1,2, Vera Fofonova 1, Sergey Danilov 1,3,4, Natalja Rakowsky 1, Sven Harig <sup>1</sup> and Karen Helen Wiltshire <sup>1</sup>**


Received: 15 April 2020; Accepted: 11 May 2020; Published: 15 May 2020

**Abstract:** A newly developed coastal model, FESOM-C, based on three-dimensional unstructured meshes and finite volume, is applied to simulate the dynamics of the southeastern North Sea. Variable horizontal resolution enables coarse meshes in the open sea with refined meshes in shallow areas including the Wadden Sea and estuaries to resolve important small-scale processes such as wetting and drying, sub-mesoscale eddies, and the dynamics of steep coastal fronts. Model results for a simulation of the period from January 2010 to December 2014 agree reasonably well with data from numerous regional autonomous observation stations with high temporal and spatial resolutions, as well as with data from FerryBoxes and glider expeditions. Analyzing numerical solution convergence on meshes of different horizontal resolutions allows us to identify areas where high mesh resolution (wetting and drying zones and shallow areas) and low mesh resolution (open boundary, open sea, and deep regions) are optimal for numerical simulations.

**Keywords:** numerical modelling; unstructured meshes; finite volume; North Sea

#### **1. Introduction**

Numerical ocean models are one of the major instruments used to understand ocean dynamics. Their area of application is usually divided into global or open ocean models (with resolutions from a few up to several tens of kilometers), regional models (that include coastal seas and whose scales are typically one to two nautical miles), and models capable of representing estuaries or certain specific processes (with horizontal scales in meters). Different basic assumptions about the physical processes to be included or excluded in a particular case are one of the reasons for these divisions; it reduces unnecessarily complicated equation-solving. For example, tides are commonly excluded from the global ocean models used in climate simulations (as in, e.g., [1]) but are dominant in the dynamics of coastal regions. Another significant reason is a limitation on horizontal discretization: models for larger domains use coarser horizontal resolutions to speed up numerical calculations, parameterizing or disregarding small-scale processes. Such limitations usually result from the finite difference method used to discretize dynamical equations in the most well-known, established models: NEMO [2], ROMS [3], MOM [4], GETM [5], and many others. While the finite difference method yields realizations quickly and easily, it is only applicable to structured meshes, making it nearly

impossible to construct meshes with variable resolution that can resolve specific areas of interest as needed (coastlines, archipelagos, and shelf breaks) but which also coarsen towards open ocean (as in, e.g., [6]). Moreover, because of computational demands, the finite difference method would not be employed to resolve big domains (ones the size of a regional sea) at the resolution needed for a coastal process (ten to hundreds of meters) unless some nesting were applied.

Studies in recent decades clearly show the necessity of combining different scales in a single model to address phenomena such as the transport of matter between the coast and the open ocean, the effect of regional processes on global ocean dynamic [7,8], or more technical questions related to regional models, such as open boundary conditions [9,10].

Nesting of two or more structured grids of different resolutions is one of the common approaches used in structured-mesh codes. A widely used one-way nesting method (meaning that information from a lower-resolution grid is transmitted to a finer grid), as in [11,12], shows good results, but the same method cannot be properly applied to explore flows in the transition zone. Two-way nesting (meaning that both grids continuously exchanging information), as in [13–15], is more complicated and incurs difficulties with smoothing and damping of signals between coarser and finer domains. This lets the most important small-scale process be filtered out or left unresolved.

The alternative to structured-mesh methods are ones designed for unstructured meshes. They are geometrically flexible and allow the resolution of the mesh to be varied (within reasonable limits), letting the mesh contour complex coastlines. Such methods are used successfully by a number of well-developed ocean models such as FVCOM [6,16,17], SCHISM [18,19], and FESOM2 [1]. Most of these models are oriented toward regional or process studies. To date, most of the experience with global large-scale application has been accumulated working with FESOM. The importance of regional and coastal model studies on various time scales is hard to overstate, be it for process studies [20,21] or climate research ([22,23]). However, in most cases, there is no dynamical link to global models [24].

With either approach, nested structured meshes or unstructured meshes, the central question is what the optimal horizontal resolution is for a specific task in the modeled region. In other words, where does the best compromise lie between the quality of the simulated dynamics and computational efficiency?

Here, we present the results of simulations performed with the newly developed FESOM-C model [25,26]. FESOM-C is a coastal branch of The Finite-volumE Sea ice–Ocean Model (FESOM2) [1]. The FESOM-C model employs hybrid unstructured meshes [27] and is based on a finite-volume discretization. It is a full three-dimensional model based on three-dimensional primitive equations for momentum, continuity, and density constituents [25]. It includes modules for the open boundary, upper boundary (interaction with the atmosphere), rivers, output, and postprocessing, all of which facilitate using this model in realistic applications. In practice, the hybrid meshes used by FESOM-C consist mostly of quadrilateral elements and include triangles only where needed to link quadrilateral cells. Compared to purely triangular meshes, this significantly increases model throughput (see, e.g., [27]). Variable horizontal resolution enables the use of coarser meshes in open sea regions but more refined ones in shallow areas to resolve important small-scale processes such as wetting and drying, sub-mesoscale eddies, or sub-mesoscale dynamics of steep coastal fronts. The general structure of the model's internal arrays, variable names, and mesh utilities is similar to that of FESOM2. Its modules for external forcing, output, and parallelization are likewise similar. However, several principal aspects render the FESOM-C model capable of representing many of the physical processes of the coastal areas, such as tides and the wetting-drying mechanism. FESOM-C also differs in that it uses a terrain-following vertical coordinate. At the same time, the closeness of both models makes establishing dynamical links between coastal and global realizations relatively easy. A detailed description of the FESOM-C model is presented in [25]. One of the goals of this study was to analyze the capabilities of the FESOM-C model as a coastal model in a realistic setup. We chose to apply the model to an area of the southeastern North Sea (the "SeNS") (see Figure 1), which is a comparatively

well-studied area. Moreover, it offers the comprehensive datasets which are essential for model setup and validation.

**Figure 1.** Bathymetry of the southern North Sea (colored contours, data from EMODnet Bathymetry portal); black lines show the mesh used for the five-year run: (**a**) a full domain; and (**b**) a zoom-in to the Cuxhaven–Helgoland area (the southeastern part of the full mesh). The background map was created with openstreetmap.org (© OpenStreetMap contributors 2019. Distributed under a Creative Commons BY-SA License).

The North Sea coast, and the SeNS region in particular, is a densely populated hub with many seaports, including Europe's biggest. Millions of livelihoods depend on the state of the North Sea. Recently-built wind farms cover a significant part of the coastal area [28,29]. The mean depth of the North Sea is about 80 m and the maximum more than 700 m, but depths in the SeNS range from 20 to 80 m. Sea-bed transformation is frequent here. The North Sea is connected to the Atlantic Ocean by the English Channel to the south and the Norwegian Sea to the north. The mean wind-driven circulation pattern is counterclockwise. Tidal dynamics are defined mainly by the semidiurnal principal lunar tide (M2). The M2 tidal wave propagates from the north and the southwest and enters the SeNS region on its western boundary. Of the three amphidromic points in the North Sea, two of them are in the SeNS. Superposition of the M2 and S2 tides causes significant spring-neap tides. The North Sea's physics are briefly described in [30]. There are several significant freshwater sources including the rivers Rhine and Elbe in the SeNS, forming strong horizontal salinity gradients. The Wadden Sea, a series of islands separated by tidal inlets and characterized by significant tidal flats and wetlands, plays an important role in the dynamics and ecological state of the SeNS [31].

As mentioned above, the North Sea, and the SeNS in particular, is a quite well-studied area. Many models have been applied to simulate its dynamics. Recent work in pre-operational modeling has focused on the southern North Sea with a combination of numerical and observational methods, as demonstrated by Stanev et al. [32]. Comparing several models, the authors found a significant difference in a nonlinear effect within the tidal dynamics of the shallow coastal zone, which they attributed in all likelihood to coarse horizontal resolutions. An evaluation [33] of two operational models with relatively high resolutions (up to 900 m) in the North Sea region, BSHcmod v4 and FOAM AMM7 NEMO, found good agreement between the models and observations of the open sea. Near the coast however, the study found a significant deviation between simulated and observed salinity. Haller et al. [33] indicated limited spatial resolution and complicated modeling in the transitional zone as weaknesses of both models. Several one-way nesting system variations, including the TRIM-NP [12] and GETM [11,34] models, have been applied successfully to the North Sea with a focus on SeNS. Gräwe et al. [11] concluded that successful modeling of tidal flow at various Wadden Sea inlets requires various minimum resolutions for some inlets. 500 m may be enough for some, but even a horizontal resolution of 200 m is insufficient for others. SCHIM, which uses finite element and finite volume discretization methods on unstructured meshes, has been successfully applied for the coupled North

Sea–Baltic Sea system with a local resolution of approximately 200 m in the SeNS [35,36]. The same model successfully applied to the Ems Estuary (Pein et al. [37]).

The main objective of this article is to demonstrate the representativeness of the latest results of the first fully realistic, three-dimensional, multi-year baroclinic hindcast simulations with the newly developed FESOM-C model. This is done by comparing model results with various observational data available for the period 2010–2015 and with those of other available models for the SeNS. Another, equally important objective is to present an application of convergence analysis of solutions for grids of different spatial resolution.

The paper is organized as follows. Section 2 describes the model setup for the SeNS. The model results for barotropic and baroclinic formulations are described in Section 3. Convergence numerical solutions for meshes with different spatial resolutions are discussed in Section 4. Section 5 offers conclusions.

#### **2. Model Setup: Southeastern North Sea**

#### *2.1. Bathymetry and Mesh*

Data from the EMODnet Bathymetry portal [38], with a resolution of about 230 m, were used to construct the model bathymetry. Scattered bathymetry data from the Wasserstraßenund Schifffahrtsverwaltung des Bundes of the German Federal Ministry of Transport and Digital Infrastructure, with resolutions of up to 1 m in riverine areas, are also employed in this region.

A low-resolution mesh (spatial resolution between 1 and 4 km with 43,318 vertices) and 21 vertical sigma layers were used to simulate the period from 2010 to 2014. This mesh was constructed by the Gmsh mesh generator [39] with the Blossom-Quad method [40] and consists mainly of quadrilaterals. Androsov et al. [25] showed that the quality of quadrilateral meshes constructed by Gmsh, even in the presence of acute angles and degenerate quadrangles, is good enough for solution convergence and the stability of FESOM-C.

#### *2.2. Boundary Conditions and Atmospheric Forcing*

For our final simulations, we used data from hydrography reconstructions based on optimal interpolation by Nú nez-Riboni and Akimova [41]. Monthly resolved data were linearly interpolated by the model on the current time step. A relaxation time parameter of 15 days (half the time of the available data resolution) in the case of propagation into the domain, and of 5 days in the case of outward propagation, were applied for temperature and salinity at the open boundary.

Sea-surface elevation at the open boundary was prescribed by amplitudes and phase for the nine (M2, S2, N2, K2, K1, O1, P1, Q1, and M4) most significant tidal harmonics in this area. Data from regional tidal solutions for the European Shelf 1/30◦ from the TPXO model ([42]) was interpolated onto open boundary locations. High grid resolution (2 km) of TPXO regional setup and improved variational data assimilation method for the shallow-water tides together with carefully selected set of tide gauges gives a good enough solution for area of open boundary. The sensitivity study of M2 wave propagation employed only M2 tidal harmonic data. These data were extracted from Danilov and Androsov [27], who modeled the full North Sea using the previous version of the FESOM-C model. Details of the open boundary implementation are provided in [25,43].

Data from the European Union's Seventh Framework Program (EU FP7) project "Uncertainties in Ensembles of Regional Re-Analyses" (UERRA) [44] constitute the surface atmospheric forcing. The atmospheric data are derived from a data assimilation method that assures its quality. The ocean model utilizes fluxes of freshwater (rain and snow), short- and long-wave radiation, surface wind, humidity, air temperature near the sea surface, and air pressure at sea level. The time resolution of the atmospheric data is 1 h; their horizontal resolution is about 11 km.

#### *2.3. Initial Conditions and Spin-Up Period*

Preliminary sensitivity studies have shown that perturbation in the initial temperature and salinity fields is compensated for after one year of simulations in such a way that the model solutions with varying initial conditions are very close to one another. Initial conditions for the one-year spin-up runs were constructed from TRIM-NP model results [12].

#### *2.4. Rivers*

Strong cross-shore salinity gradients in shelf areas like the SeNS are mainly determined by the supply of fresh water from rivers. Temperature and observed daily river-runoff, put together by Kerimoglu et al. [45] from Radach and Pätsch [46], were used to prescribe freshwater supply. A salinity of 0.1 [psu] was used for river water to avoid numerical instabilities due to frontal dynamics near fresh water sources. In total, nine freshwater sources were prescribed (see Table 1).


**Table 1.** Freshwater sources in the current setup based on daily observed data.

During the simulation period, several "flood" events comprising a significant increase in water discharge were observed during the winters of 2011, 2012, and 2013 and the summer of 2013. Details on the effects of the 2013 summer flood event are described in [47].

#### **3. Simulation Results**

In this section, we provide a basic validation of the simulation carried out on a mesh with variable resolutions of about –4 km in the SeNS area for the five years from 2010 to 2014. As previously mentioned, the SeNS is the most heavily observed area, with the amount of data increasing significantly over the last few years as new instrumentation has been developed. Several databases, such as EMODnet and COSYNYA, contain a significant part of this data. Nevertheless, a consistent database and a method for model validation in this area are lacking.

#### *3.1. Tidal Dynamics*

Tides are one of the main driving forces in this area. Accurately representing them is one of the most critical factors in any successful description of the coastal dynamics.

Previous work has demonstrated that the model can accurately reproduce tidal dynamics ([25–27]). To test how the current model setup performs at reproducing the main tidal harmonics in long-term simulations, we analyzed observed sea-level height at several regional stations by extracting the amplitudes and phase of nine harmonics for comparison with model results. We also performed a sensitivity run with only the M2 harmonic prescribed at the open boundary.

#### *3.2. M2 Tide*

To test the ability of the model to reproduce the main tidal wave (M2 harmonic) in the SeNS domain, we constructed an additional 2D barotropic setup (by switching off the baroclinic part) with only the elevation from the M2 tide wave prescribed at the open boundary. To get a clear tidal wave and make the analysis simpler and more accurate, we disregarded atmospheric forcing for this 2D setup. Analyzed upon onset of an equilibrium regime after several days of simulation, the results of this simulation were compared with observations, as shown in Figure 2. The M2 tidal constituent propagates along the coast of the North Sea as a Kelvin-type wave. It enters the model domain at the western boundary and propagates eastward to the Elbe River estuary and then northward along the coast. The SeNS area is characterized by two amphidromic points (points of zero amplitude): one in the southwest around 3.5◦ E, 52.5◦ N, and the other in the north around 5.5◦ E, 55.2◦ N. Both of the M2 wave's amphidromic points are well represented in the model (Figure 2). The amplitudes and phases were compared to observed values provided by Ole Baltazar Andersen (personal communication, 2008). The amplitude and phase accuracy is characterized by the total vector error ([25]):

$$\mu = \frac{1}{N} \sum\_{n=1}^{N} \left( \left( A\_{\ast} \cos(\phi\_{\ast}) - A \cos(\phi) \right)^{2} + \left( A\_{\ast} \sin(\phi) - A \sin(\phi) \right)^{2} \right)^{1/2},\tag{1}$$

where *A*∗, *φ*<sup>∗</sup> and *A*, *φ* are the observed and simulated amplitudes and phases, respectively, at *N* stations. For the current model setup and 53 observational stations, the total vector error is 0.21 m, compared with a maximum of wave height of 2 m. This can be regarded as a good result. The most significant error is simulated near the coast. Most of the discrepancy can be explained by uncertainty in the model's bathymetry and bottom-drag parameterization. In general, the model reproduces the observed values reasonably well but with some exceptions. Stations in the area of Cuxhaven (to the west of the domain) show a smaller amplitude than the observational dataset provided by Ole Baltazar Andersen. However, in the experiment with nine prescribed tidal harmonics at the open boundary, the stations in this area (Helgoland and Cuxhaven) are reproduced well by the model, even the M2 amplitudes and phases (Figure 2).

#### *3.3. Main Tidal Harmonics on Long Time Scales*

The effect of neap-spring tide variability is known to be important for coastal dynamics. A number of previous studies, both observational and model-based, have shown the importance of tidal harmonics besides M2; in the studied area in particular, tidal non-linearity also plays a role in the dynamics of the coastal areas ([32,48,49]). We prescribed nine tidal harmonics at open boundaries for our 3D baroclinic run for the years 2010–2014. We compared the amplitudes and phases of the nine main ones simulated with FESOM-C, and as observed at four stations, with results from the tidal solutions of the TPXO model, European Shelf 1/30◦ version. The stations' positions are indicated by the black dots in Figure 2a,b. The observed amplitude and phase values are based on the years 2010–2014 and calculated using the uTide python module ([50]). The stations used to validate the model were selected for their location at points of interest in the study of tidal dynamics and for providing observed values that cover the simulated period. Two stations, K13a3 and Hoek van Holland, lie near one of the amphidromic points of the M2 tidal wave. The Helgoland and Cuxhaven stations are located near the position where the M2-only simulations yielded the maximum deviation from observations. The Hoek van Holland and Cuxhaven stations are both land stations. The amplitudes and phases simulated by FESOM-C correspond well to the observations; comparing the M4 tidal constituents, we find them just as good as or better than the TPXO model. The Hoek van Holland station is exceptional: here, FESOM-C shows higher amplitudes for M2 and M4. Its vicinity is poorly resolved in the current version of the model setup, with mesh vertices approximately 4 km apart, so bathymetric errors play the most significant role here. Moreover, data from this area are assimilated into the TPXO model. However, the results at the Cuxhaven station, where FESOM-C has significantly higher resolution than TPXO, show the opposite: FESOM-C reproduces all the tidal harmonics very well. Direct comparison of FESOM-C results with observations shows high correlations between modeled and observed values. The standard deviations (STD) are quite close for all stations (see Table 2).

**Figure 2.** Comparison of the modeled and observed characteristics of the M2 tidal wave for: the amplitude, in meters (**a**); and the phase, in degrees (**b**). The color maps present model results; colored circles correspond to observations. (**c**,**d**) Scatter plots for amplitude and phase, respectively, for the entire domain. The numbers in the panels are the ID numbers of the stations. The total vector error is 0.21 m. The black circles and text in (**a**,**b**) indicate the positions and names of the stations in Figure 3.

**Figure 3.** Amplitudes (**a**,**c**,**e**,**g**) (meters) and phases (**b**,**d**,**f**,**h**) (degrees) of the nine main tidal harmonics at four stations in the modeled area. Green shows the calculated ([50]) values based on observed sea surface elevation during the years 2010–2014. Blue shows the regional (European Shelf 1/30◦) tidal solutions of the TPXO tidal model. Cyan shows FESOM-C model results processed in a way similar to the observations.


**Table 2.** Correlation between observed and simulated sea surface height (ssh) at four stations. Standard deviation (STD) in time arrays of observed and modeled values of ssh.

#### *3.4. Surface Salinity and Temperature over 2010–2014*

Modeled salinity and temperature are verified on the basis of data from ICES (Figure 4), COSYNA (maintained by the Helmholtz Zentrum Geesthacht and accessible at www.cosyna.de ([51]), and EMODnet (Figures 5–7) databases. The model captures the observed ([52]) lateral salinity gradient (Figure 4a) reasonably well. Due to residual barotropic currents, a strong horizontal salinity gradient forms along the coast. Figure 4b,c, respectively, provide a comparison of modeled surface salinity and temperature with observed data from the ICES database for the years 2010–2014. With some exceptions, the observations are from locations outside areas with a strong horizontal salinity gradient, in the open area of the simulated domain, where the vertical structure of the water masses is more susceptible to variability associated with a seasonal thermocline. The corresponding Pearson correlation coefficient (Cor.) and root mean square difference (RMSd) are indicated at the top of the plots. Simulated surface salinity and temperature correlate well with observations. The high-temperature correlation of 0.99 can be explained by the seasonal cycle. The relatively small RMSd for both temperature and salinity attests to the model's ability to reproduce the overall dynamics and thermohaline structure of the simulated region.

**Figure 4.** (**a**) Mean sea surface salinity for the years 2010–2014. The black dots indicate positions of data from the ICES database; black-with-red circles indicate positions of stations giving temperature and salinity time series with the corresponding names. (**b**,**c**) Sea surface salinity and temperature from FESOM-C (y-axis) compared one-to-one with values from the ICES database (x-axis) with corresponding correlation coefficients (Cor.) and root mean square difference (RMSd).

#### *3.5. Time Series of Temperature and Salinity*

The southern North Sea area is rich in observational data, making this area interesting for model calibration and validation. We validated the model using an automated validation system introduced in the new model version. A special model-structure output module provides output at discrete predefined stations for direct comparison to the observed time-series collected in various databases. The collection is much larger than the set of stations dealt with above; we only show comparisons for stations with interesting dynamics from which continuous measurements are available.

In Figures 5–7, a time-series of simulated temperature and salinity for the years 2010–2015 is compared with observational data from several autonomous stations. Most of the data taken from the Emodnet and COSYNA databases underwent automatic quality control. However, some data remain doubtful, for example a somewhat high temperature reading at the Sylt station during the winter of 2012–2013 (Figure 5c). Filtering and cleaning observations is beyond the scope of this paper. Nevertheless, most of the data are trustworthy and useful for direct comparison.

**Figure 5.** Observed (grey dots) and modeled (blue line) sea surface temperature at three stations: (**a**) Den Helder; (**b**) K13a3; and (**c**) Sylt. Station positions are indicated in Figure 4.

**Figure 6.** Comparison of observed (grey dots) and simulated (blue line) temperature at two stations: (**a**,**c**) the Helgoland station at 1 and 10 m depth, respectively; and (**b**,**d**) the UFS Deutsche Bucht station at 20 and 30 m depth, respectively.

The seasonal cycle of surface temperatures at three stations, namely K13, Den Helder, and Sylt (Figure 5), as well as at Helgoland (Figure 6a), was captured well by the model. The Den Helder station (Figure 5a) is situated on the first Wadden Sea inlet and experiences intense water exchange between the inner Wadden Sea and the open sea; in part, this explains the high temporal variability in observed and simulated temperature. The K13a3 station (Figure 5b) is situated close to the model domain's open boundary (here representing the open sea) and is less affected by coastal processes. The Sylt station (Figure 5c) is situated close to one of northern Wadden Sea bays used for various model test cases and near a salinity gradient under the influence of freshwater. The Helgoland station is characterized not only by highly variable salinity and a strong lateral salinity gradient but also by a significant bathymetric gradient. In general, the model captured observed surface temperature

dynamics, including cyclical-seasonal and local effects, very well at these stations. The model generates a lower autumnal temperature at the K13a3, but there is no such effect at the other stations.

**Figure 7.** Comparison of observed (grey dots) and simulated (blue line) salinity at three stations: (**a**) data from the Helgoland station at 1 m depth; (**c**) data from the Brouwershavensegat station at 1 m depth; and (**b**,**d**) data from the UFS Deutsche Bucht station at 6 and 30 m depth, respectively.

Temperature validation for deeper layers is shown in Figure 6 for the two stations Helgoland (1 and 10 m depth) and UFS Deutsche Bucht (20 and 30 m depth). UFS Deutsche Bucht is not far from Helgoland; however, there is a deep trench leading to the open sea from the latter. Here, temperature undergoes a pronounced seasonal cycle with less seasonal variability than at the surface. The model exhibits greater deviation from observations in deep layers than at the surface. It represents observations reasonably well in general, albeit with local differences up to several degrees.

Figure 6 represents observed and simulated salinity time-series at three stations: Helgoland (1 m depth), UFS Deutsche Bucht (6 and 30 m depth), and Brouwershavensegat. The Brouwershavensegat station (1 m depth) is located near the open boundary and is affected by water from the Rhine. There are no pronounced seasonal salinity dynamics in the SeNS area ([11]). Salinity changes at these stations are determined mainly by changes in wind-driven currents and freshwater supply from rivers. Observational data variability on small time scales (of days) is generally higher than in model results; this is to be expected under the current setup, with 21 sigma layers and horizontal resolutions of up to 1 km that will not resolve river plumes and freshwater lenses properly. While the model does not track rapid changes in observed salinity, it does capture common dynamics well. Drops in salinity due to extreme flooding events such as the one in 2013 are clearly seen in observation and reproduced by the model, as shown in Figure 7a–c. The model shows a drop salinity at 30 m depth during 2013, unlike observation (see Figure 7d). This could be related to the current setup's rough vertical resolution and parameterization of vertical mixing. The most significant difference between observation and the model appears for the year 2011.

#### *3.6. Salinity and Temperature, Ferry Lines*

The near-shore area of the SeNS is characterized by strong lateral density gradients defined by the salinity gradient, mainly due to fresh water supply from rivers ([11]), and by the temperature gradient under the effect of the different temperature dynamics of shallow and deep areas. Such gradients may play a significant role in near-coastal dynamics ([34,53,54]). Temperature and salinity data collected

by the Operational Systems department of the Institute of Coastal Research Helmholtz-Zentrum Geesthacht, by way of FerryBox systems installed on several ferries connecting Cuxhaven–Helgoland (CH), Büsum–Helgoland (BH), and Cuxhaven–Immingham (CI) [47,55,56], offers a way to verify the model both in near-coastal areas (BH and CH) as well as farther offshore (CI). These ships' areas of operation are indicated by the three ellipses in Figure 8, superimposed over their respective paths indicating salinity values. Raw data for analysis were taken from the COSYNA database. The CH and BH ferries both call at the island of Helgoland, where mean salinity is about 34 PSU. The Port of Cuxhaven (serving CH and CI) is situated at the western side of the mouth of the river Elbe, and the port at Büsum is about 30 km north of the Elbe estuary. Both are located in the region of freshwater influence from the Elbe, with horizontal salinity gradients up to 0.45 PSU/km and tidal amplitudes in excess of 1.5 m, with a broad wetting-and-drying area. We used data from the CI ferry between 3◦ E longitude and Cuxhaven; the rest lay outside the model domain. One axis in the figures tracks longitude; the sailing routes run mainly east–west. The CI ferry altered its route over time, as is reflected in the mean values (Figure 9).

**Figure 8.** Sea-surface salinity measured by three ferry lines across their paths used for comparison in Figures 10–12. Greyscale background shows bathymetry. The color scatter plots show observed salinity for the period 2010–2014. The three ferries' areas of operation are indicated by ellipses: (**b**) red for Cuxhaven–Helgoland (CH), blue for Büsum–Helgoland (BH), and (**a**) green for Cuxhaven–Immingham (CI).

**Figure 9.** Mean sea-surface salinity from the three ferry lines (black lines) and corresponding model values (red lines). Mean values are shown by solid lines; mean values ± one standard deviation are shown by dashed lines. Ferries: (**a**) Cuxhaven–Helgoland; (**b**) Büsum–Helgoland; and (**c**) Cuxhaven–Immingham.

The data sampling rate varies from 10 to 50 s and depends on time and route; the distance between data samples thus depends on a ship's speed and the sampling rate. The distance between neighboring measurements across the dataset used here varies from 80 to 400 m. Both the spatial (80 m) and temporal (10 s) resolutions are much higher than in the model (1000 m in coastal areas and 60 s, respectively) and especially higher than in the three-dimensional model output (with the same spatial resolution but 1.5-hourly mean values). The number of model-3D-output snapshots is sharply limited by the space required to store the output data and the performance of long-term memory (HDD and SSD), and could not be significantly changed. Observed data are also scattered across time and space, which limits the ability to save model output at exact times and positions as doing so would mean storing more than 2 million scattered data points that are not yet realized in the model. Subsequent versions of the model will include this feature. The distinction between the sampling rate and model output introduces a discrepancy into any direct comparison. A temperature difference of 0.2–0.4 ◦C arises due to the time shift. Comparing coarse horizontal model resolution to observation leads to a situation in which several observational points correspond to a single model point. Due to salinity gradients, and accounting for tidal currents of up to 2 m/s, errors in salinity can approach ±2 PSU. The more transects are sailed, the less the mean values will differ due to random errors such as in spatial and temporal resolution; however, errors in deviation will remain. Greater differences are expected in near-coastal areas because the water masses there vary much more. Salinity goes from up to 35 PSU near Helgoland to brackish near land, with salinities between 5 and 10 PSU. Direct comparison among different instruments (FerryBox, water sample analyses, OSTIA satellite data, and MARNET stations) were provided by Haller et al. [33], Grayek et al. [57]. Haller and colleagues [33] reported an error of 0.79 RMS in salinity between FerryBox data and laboratory water sample analyses. Grayek and colleagues [57] found that temperature differences between FerryBox and MARNET station readings can be the result of heating of the measured water in the ferry.

The difference in the density of water masses is known to determine the dynamics of the baroclinic system. Similar to Gräwe et al. [11], in Figures 10 and 11, we give temperature and salinity anomalies for the three ferry lines. Apparently, both density constituents (for salinity, see Figure 10; and, for temperature, see Figure 11) play a significant role in the formation of the density gradient (Figure 12). This comparison allows us to evaluate the potential source of error in the model. Anomalies were calculated separately for each ferry section, eliminating seasonal differences. However, comparing only the anomalies in the model does not inform us about absolute values. The statistics for this comparison are shown in Table 3 and Figures 9 and 13 for absolute values and anomalies.

**Figure 10.** Salinity anomalies from the three ferry lines and corresponding model results. The observational data are on the left and the corresponding model results are on the right of each panel. Ferries: (**a**) Cuxhaven–Helgoland; (**b**) Büsum–Helgoland; and (**c**) Cuxhaven–Immingham.

**Figure 11.** Temperature anomalies from the three ferry lines and corresponding model results. The observational data are on the left and the corresponding model results are on the right of each panel. Ferries: (**a**) Cuxhaven–Helgoland; (**b**) Büsum–Helgoland; and (**c**) Cuxhaven–Immingham.

**Figure 12.** Density anomalies from the three ferry lines and corresponding model results. The observational data are on the left and the corresponding model results are on the right of each panel. Ferries: (**a**) Cuxhaven–Helgoland; (**b**) Büsum–Helgoland; and (**c**) Cuxhaven–Immingham.

**Figure 13.** Distribution of salinity in measurements from the three ferry lines (blue) and the model (green). Ferries: (**a**) Cuxhaven–Helgoland; (**b**) Büsum–Helgoland; and (**c**) Cuxhaven–Immingham.

**Table 3.** Comparison of salinity (S), temperature (T), and density (*ρ*) between FerryBox and FESOM-C. RMSE, Root Mean Square Error of anomalies. The numbers in brackets are based on absolute values.


Basic statistics for the available data and comparison with the model are shown in Table 3. The total number of measured points from over 2000 ferry transects is more than 2,000,000. Data from the CI and BH ferry lines cover the whole modeled period of 2010–2015; the gaps (the white area in the figures) are seasonal, reflecting winter (BH) or summer (CH). The CI ferry line has fewer gaps for 2010–2012. The second part of the simulated period is not covered by data from the CI ferry except for some data from the end of 2014. The CI ferry line traces a much longer transect compared to the other two and mostly represents water masses with a high salinity of more than 34 PSU (see Figure 13, blue bars). The BH and CH ferry lines mostly sail near the fresh water influence area of the river Elbe, where salinity varies 10–35 PSU for CH and 20–35 PSU for BH. The RMS error of anomalies is smaller in general than the corresponding RMSE of absolute values shown in brackets. Gräw and colleagues [11] compared results of the well-established GETM model with BH ferry data for the 2009–2011 period. The horizontal resolution of the GETM grid presented by Gräwe et al. [11] is 200 m, that is, five times finer than the present setup. Gräwe et al. [11] reported RMSEs for salinity of 1.15 PSU, 0.64 ◦C for temperature, and 1.32 kg/m<sup>3</sup> for density. Haller and colleagues [33] compared the results of two three-dimensional hydrodynamic models, BSHcmod and AMM7, with the CI ferry line the years 2009–2012. BSHcmod v4 is a three-dimensional operational model with a two-way nesting approach and a finest resolution of 900 m. AMM7 is a one-way nesting, operational model based on the NEMO model that includes assimilation of in-situ observations at a 7-km horizontal resolution. Haller et al. [33] calculated salinity RMSEs of 0.68 and 1.1 PSU and temperature RMSEs of 0.68 and 0.44 ◦C for the BSHcmod v4 and AMM7 models, respectively.

The results of FESOM-C are within common statistics and clearly show better agreement in the open sea. However, comparison with the BH ferry data indicates a problematic area.

The anomalies of temperature, salinity, and density are shown in Figures 10–12, respectively, for the three ferry lines CH, BH and CI. The seasonal cycle is easily visible in the shift in positive and negative temperature anomalies from summer to winter. Such a shift (both in time and position) is reproduced well by the FESOM-C model on a small spatial scale of about 60 km (Figure 11b) as well as on longer spatial scales of 300 km (Figure 11c). Some of the small features, such as surface warming at 8◦ E during the autumn of 2014 (see Figure 11c), are also captured reasonably well by the model. The RMS difference for the CI ferry line is close to the comparison made with data from the ICES database (Figure 4). The difference between the mean observed and the modeled temperature varied from 0.4 to 0.8 ◦C and increased with rising temperatures. The temperature difference increases significantly towards the coast in shallow waters. Several possible reasons for this mismatch are worth mentioning; although their importance is known, the current setup still does not implement them. The depth to which short-wave solar radiation penetrates varies significantly in this region and depends heavily on suspended matter. There is a well-known steep gradient in suspended matter along this region's coast. The lack of feedback with the atmospheric model used at the upper-boundary ocean model is an issue too, since the atmospheric model assumes surface-water temperatures different from those simulated by FESOM-C. That leads, in turn, to incorrect long-wave radiation ([58]). Inaccuracy in sea-bed albedo, together with wave effects, might also modify surface heat flux. These shortcomings of boundary condition formulation for tracer equations will be taken into account in the next version of the model.

Unlike temperature, salinity, and density do not show a pronounced seasonal cycle. However, salinity exhibits steep, more pronounced offshore gradients. The model reproduces salinity and density anomalies reasonably well, as it does with temperature. Most of the differences appear in shallow areas with maximum lateral gradients. Variability in salinity is generally close to what is observed, except for one area near Helgoland (the western port of the CH and BH ferries; see Figure 9a,b). Here, the model's standard deviation is about two times smaller than observations while mean salinity is close to observations. The temporal dynamics are captured reasonably well by the model at the Helgoland station (see the salinity time-series in Figure 7a). Simulated mean salinity starts to deviate significantly from observation towards the coast in the area north of the mouth of the Elbe River (the BH ferry line). The model reproduces observed salinity dynamics west of the Elbe (the CH ferry line) well.

Several factors may explain the model's deviation in the area of the BH ferry line. The tidal wave propagates from west to east up to the Elbe River and then northward along the coast. The tidal dynamics near Cuxhaven (the eastern port of the CH ferry) is reproduced well by the model. However, north of the Elbe River, the model shows a significantly lower amplitude for the M2 tidal wave (see Figure 2, Station 12). Stations 14 and 15 in the same figure, lying further north, do not show such a significant deviation. Areas of wetting and drying play a significant role in this region's dynamics. A significant near-coastal area falls dry during low tide. Errors in the model's bathymetry—the area north of the Elbe River is deeper—and coarse resolution in the wetting-and-drying areas, together with uncertainties in bottom drag, lead to higher Elbe-water transport along the coast. However, most likely, it should rather propagate in a north-westerly direction, towards Helgoland. A sensitivity study with higher resolution near the coast and a more accurate representation of the bathymetry improved both barotropic and baroclinic model dynamics in this area. The results of these studies are partly presented in Section 4.

#### *3.7. Vertical Structure, Gliders*

For the 2010–2014 period, the Operational Systems department of the Institute of Coastal Research Helmholtz-Zentrum Geesthacht performed several surveys, with two gliders called "Sebastian" and "Amadeus", in the area of Helgoland and towards the open sea. Data from these surveys (which cover an area not reflected in the FerryBox data; see Section 3.6) are available in the COSYNA database. With high vertical and special resolutions, these data offer a unique opportunity to validate the model. Modeled salinity and temperature are compared with glider data in Figures 14 and 15, respectively, in which bathymetry with glider path is shown at right, model data at top left, and glider data at bottom left of the panels showing the two surveys. Time of measurement is shown on the map by color, from red (beginning of survey) to yellow (end of the survey); this corresponds to the time axis on the plots. Every glider measurement is shown by a colored circle. The model output was interpolated to match the glider data in time and space by the "nearest neighbor" method. For every glider data point, one profile corresponding to the time of measurement is extracted from the model output, to exclude duplicating profiles. The "nearest neighbor" method introduces some inaccuracy and moreover allows only point-by-point comparison between the model and observations. The measured data points total about 60,000 (winter transect) and 26,000 (summer transect), vastly fewer than the FerryBox data, but are distributed in 3D space and so require three-dimensional interpolation. 3D linear interpolation of the model results on an unstructured mesh, as was done for the 2D FerryBox data, was dispensed with as time and memory consuming. Furthermore, the 3D model output that formed the basis for comparison was saved at 1.5-h intervals. The time resolution of the glider data is about 1 min, so the horizontal resolution of the observed data is thus significantly higher than in the model. The 1.5-h time resolution of the model output introduces an additional discrepancy as against observation due to tidal phase shifts. The water masses may have shifted about 5 km horizontally, given a current of 1 m/s. These uncertainties in data comparison should be taken into account during future analysis. However, the resulting figures are nevertheless useful to point out the advantages and disadvantages of the model.

The results of two glider transects from the winter of 2011 and the summer of 2013 are shown in Figure 14 (salinity) and Figure 15 (temperature). The 2011 winter glider path was from near Helgoland (9 February 2011) northwest for about 60 km and back (22 February, 2011). Quality checked salinity data are available only for the second part of the transect, starting on 15 February, 2011. Both the temperature and the salinity figures show a well-mixed water column in observations and model results. Shortly before and during the winter transect, winds were strong as compared to the summer transect. Mean wind speed was about 8 m/s (and up to 16 m/s) with mean air temperatures around 0 ◦C near Helgoland. This explains the mixed water column and slight overcooling at the surface in the model. Horizontal temperature (lower temperature near Helgoland) and salinity (increase salinity towards the open sea) gradients are reproduced reasonably well by the model. Observations show wavy variability in both the temperature and salinity profiles, with a period similar to the M2 tidal wave, which is probably the effect of a tidal dynamic. The model captures a similar dynamic well. Unlike the winter situation, the summer transect (see Figures 14b and 15b) shows strong vertical stratification in the temperature fields. Mean wind speed on the summer transect was 5.5 m/s (and up to 12 m/s), mean air temperature about 18 ◦C, and short-wave solar radiation was about four times higher than in winter, determining the extant strong thermocline. While the model shows similar sea-surface and near-bottom temperature dynamics, it does not capture the sharp vertical gradient. The simulated temperature and salinity profiles are smoother than the observed ones. The model does not resolve surface freshwater plumes near Helgoland. Vertical gradients in modeled salinity profiles are much less pronounced, and smoothed, compared to the vertical temperature structure. The current setup has 21 sigma layers, which may not be enough to capture steep vertical gradients. A sensitivity test with more layers shows some improvements in model results. Improvements in vertical turbulence schemes are also required.

**Figure 14.** Comparison of model salinity with observed data from gliders: (**a**) February 2011; (**b**) August 2013. Colored circles on bathymetry maps (right) show glider paths; color indicates position times from red (beginning of the survey) to yellow (end of the survey). The filled contours at upper left are model results (21 sigma levels, 2010–2014 run). The scatter data at bottom left are from gliders (COSYNA database ([51])). The glider data time resolution is about 1 min, or 20 m horizontal resolution. The spatial resolution of the model is between 1 and 4 km.

**Figure 15.** Comparison of model temperature with observed data from gliders: (**a**) February 2011; (**b**) August 2013. Colored circles on the bathymetry maps (at right) show glider paths; colors indicate position times from red (beginning of the survey) to yellow (end of the survey). The filled contours at upper left are model results (21 sigma levels, 2010–2014 run). The scatter data at bottom left are from gliders (COSYNA database [51]). The glider data time resolution is about 1 min, or 20 m horizontal resolution. The spatial resolution of the model is between 1 and 4 km.

#### **4. Discussion, Mesh Resolution**

Comparing the model results with observations, we find the largest discrepancy in shallow areas near the coast, especially in wetting-and-drying areas. Previous work, as well as sensitivity studies with the current model, have shown that refining the horizontal resolution of the mesh significantly improves model results. At the same time, it is not always obvious what horizontal resolution different regions require. To find the optimal resolution, we performed several sensitivity simulations.

#### *Solution Convergence on Different Meshes*

One of the important stages during preparation is selecting the optimal mesh resolution for the modeled region. We understand optimal mesh resolution as a compromise between computational efficiency and the quality of the simulated dynamics. A preliminary calculation on the sequence of meshes allows us to estimate the convergence of numerical solutions. We constructed three different meshes for our experiments. All were generated by the Gmsh mesh generator ([39]). The first one (m8) has a spatial resolution of between 4 and 1 km and 43,318 vertices. The resolution of the second mesh (m5) varies from 2.2 to 550 m and has 134,858 vertices. The third mesh (m3) has minimum and maximum cell sizes of 250 m and 1.6 km, respectively, and 235,283 vertices. All meshes have 21 non-uniform sigma layers in the vertical direction (refined near the surface and bottom). The wetting/drying option was turned on. To test the code's sensitivity to the mesh resolution, we computed barotropic, tidally driven circulation in the SeNS for two atmospheric scenarios: one without a strong wind effect for the full tidal period ("weak wind"), and the other with a "strong wind" component. Discrete sea-surface height (SSH) values and two components of velocity (*u*, *v*) were accumulated in the mesh vertices, and then these values were linearly interpolated onto the vertices of the coarse mesh (m8). Next, we determined the values of the different *δζ*, *δu*, and *δv* among solutions on the meshes. We performed the comparison for the full M2 tidal cycle, as shown in Figure 16, which presents the histograms of the differences.

For the solutions on the m8 and m5 meshes for the elevation values (the tidal wave maximum exceeds 3.5 m), only 30.6% of the points agree within a range of ±1 cm for the "strong wind" experiment; 31.7% agree in the "weak wind" scenario (Figure 16(top)a,b). For the ±2 cm interval, the proportion of points with the same solution more than doubles, while exceeding 67% for both scenarios. A similar situation occurs when we compare horizontal velocity components on these two meshes (Figure 16(bottom)a,b). In the range of ±1 cm/s, the number of identical points for the

*u*-component is approximately 43% in both scenarios but slightly exceeds 48% for the *v*-component (the maximum horizontal velocity is about 280 cm/s). For a ±2 cm/s interval, the number of identical solutions is about 80% for both velocity components.

Comparing the more detailed m5 and m3 meshes, we see a significant improvement in convergence of the results, with the exception of the sea-surface height comparison in a "weak wind" scenario (Figure 16(top)c,d). Here, the number of points falling within the ±1 cm limit is only 27.9%; in a "strong wind" experiment, the number of such matches is almost twice as large (55.6%). Increasing the match interval up to ±2 cm, the number of matches increases significantly, to 73.3% for a "weak wind" and 85.1% for a "strong wind" scenario. A significant improvement in convergence of the results is also apparent in the horizontal velocity field (Figure 16(bottom)c,d). Thus, in the range of ±1 cm/s, the number of matches in the "weak wind" experiment is 55% for the transverse velocity component and 62% for the longitudinal velocity component. In the "strong wind" case, the convergence of results is greatly improved, exceeding 70% for both components of the velocity vector. In the range of ±2 cm/s, the convergence of solutions tops 90% for the velocity components of the two experiments.

As is apparent from the foregoing results (see Figure 16), the experimental "strong wind" simulations negate the poor spatial resolution somewhat, especially in the shallow-water zone where mesh resolution errors are highest. Either onshore winds increase the thickness of the water, or offshore winds drain the shallow-water zones faster; the result is smoother errors on meshes of different resolutions.

**Figure 16.** Histograms of the difference between solutions at different mesh resolutions: (top) sea-surface elevation; (bottom left) the *u*-component of velocity; and (bottom right) the *v*-component of velocity. (**a**) Meshes m8 and m5 ("weak wind"); (**b**) Meshes m8 and m5 ("strong wind"); (**c**) Meshes m5 and m3 ("weak wind"); and (**d**) Meshes m5 and m3 ("strong wind").

As can be seen from the results of the spatial comparison (Figure 17), the maximum difference falls in zones of drying and of minimum depth. In detailing the coastal zone, the wetting-and-drying processes differ, to an extent, from the solution on coarser meshes. The difference between the solutions on the coarse and middle-resolution meshes in the "weak wind" case reveals a certain difference in the solutions for the deep-water part of the region precisely in the amphidromic zone of the M2 wave (Figure 17a). The difference never exceeds 0.2 cm, meaning there is a slight shift of the amphidromic point for the solutions on the two meshes. The solutions on the intermediate and on the most detailed meshes do not show this difference in the amphidromic zone (Figure 17b).

A different situation arises in the "strong wind" scenario: Two zones of maximum difference appear in a deep-water part of the SeNS. Their locations stem from the region's bathymetric features, namely a small underwater sill. At the same time, a coarse-mesh solution produces somewhat underestimated results in elevation at localized areas not exceeding 0.25 cm (Figure 17c). There is practically no difference between the solutions on the intermediate and the most detailed meshes (see Figure 17d). This analysis shows that, for model simulations, the m5 mesh will yield solutions of optimal quality but also that a coarser mesh (m8) solution will not introduce significant errors in the model results.

**Figure 17.** The spatial difference for SSH between solutions on different meshes with different atmospheric scenarios: (**a**) at *δζm*<sup>8</sup> <sup>−</sup> *δζ<sup>m</sup>*5; (**b**) at *δζm*<sup>5</sup> <sup>−</sup> *δζm*<sup>3</sup> in "weak wind"; (**c**) at *δζm*<sup>8</sup> <sup>−</sup> *δζ<sup>m</sup>*5; and (**d**) at *δζm*<sup>5</sup> <sup>−</sup> *δζm*<sup>3</sup> in "strong wind". Red ellipses denote zones of maximum differences in the deep-sea part of the North Sea.

#### **5. Conclusions**

The first fully realistic, three-dimensional, multi-year baroclinic hindcast simulations were performed with the newly developed FESOM-C model and were comprehensively validated in the area of the southeastern North Sea. The FESOM-C model, developed mainly for the coastal region, employs mixed unstructured-mesh methods ([25,27]) and a finite-volume discretization. FESOM-C is a fully resolved three-dimensional model based on primitive equations for momentum, continuity, and density constituents ([25]). Well-developed modules for the open boundary, the upper boundary (interaction with the atmosphere), rivers, output, and postprocessing allow this model to be used

for realistic simulations. The present work uses hybrid meshes that combine quadrilaterals and a small number of triangles. Such meshes support zooming-in to an area of interest (in this case, for study of the Wadden Sea and its estuaries) even as mesh resolution towards the open sea coarsens significantly. This variable horizontal resolution enables a more efficient use of computational resources while refinement in shallow areas resolves important small-scale process (such as wetting and drying, sub-mesoscale eddies, and the dynamics of steep coastal fronts). Using meshes composed mostly of quadrilateral cells allows a significant increase in calculation rates. Proper representation of physical processes, both near shore and in the open sea, makes it possible for the model to represent their dynamics.

Overall, our validation shows that the FESOM-C model reproduces the physical dynamics of the southeastern North Sea reasonably well.

The FESOM-C model reproduces the tidal dynamics, one of the most important dynamic components of the area in question, well (on tidal dynamics, see Section 3.1). It captures the amplitudes and phases of the main tidal harmonics well compared to other solutions, and modeled sea surface height for the years 2010–2015 agrees well with observations. The model also reproduces mean and seasonal horizontal and temporal distributions of temperature and salinity reasonably well (see Sections 3.4 and 3.5). An analysis of temperature, salinity, and density readings from three FerryBoxes shows that the model can reproduce the water mass characteristics of the open sea very well, with somewhat large (but explicable) inaccuracies in the coastal zones.

An analysis of a comparison between glider data and simulated three-dimensional temperature and salinity yields estimates that allow us to evaluate the model's skill at modeling the southeastern North Sea. The vertical temperature and salinity distribution indicates that the vertical turbulence scheme needs improvement.

New developments in model output and postprocessing methods allow model results to be validated against observations on unstructured mixed meshes. Detailed validation includes various observational datasets from different autonomous instruments, such as FerryBoxes, gliders, and buoys, which are spread across time and space. The present work uses data with much higher local spatial and temporal resolutions (up to several seconds and meters) than most models, illustrating high natural variability in the coastal area.

The fact that horizontal salinity and temperature gradients as well as frontal dynamics are so well captured demonstrates that the cell-vertex (finite volumes) discretization method, used here with hybrid meshes, is capable of realistic application in this region, where the dynamics of steep gradients are a crucial and quite often difficult modeling issue. The present work should allay skepticism about unstructured-mesh coastal models as being "too dissipative" in character. In fact, they are on par with the more traditional models formulated on structured rectangular grids. The combination of coarse resolution towards the open boundary with a sharper focus on the coastal area of the Wadden Sea region offers a unique opportunity to simulate a region of interest at relatively high resolution with significantly reduced computational cost. The final computation was performed with only 24 CPUs using OpenMP parallelization. Such resources nowadays are comparable to state-of-the-art laptops and smartphones and make simulations of complicated 3D realistic cases possible without operational difficulties and expensive supercomputers. The relatively small (in terms of horizontal size and resolution) mesh used in the current setup could be extended to a larger area without significant impediments. Increasing mesh resolution by a factor of 4 leads to the number of nodes increased by a factor of 16 for a regular rectangular grid. This changes depending on the shape of the calculated domain. There is no difference between the structured and unstructured meshes in the case of the open ocean. In the proposed work, similar increase in the resolution from robust mesh (m8) to fine mesh (m3) leads to an increase in the number of nodes by a factor of 5.5. The scalability of the FESOM2 model was recently studied by Koldunov et al. [59]. In [59], the authors demonstrated that the FESOM2 model is competitive tools for a high-resolution climate modeling. In turn, FESOM-C uses the mixed (predominantly quadrilateral) meshes. Approach based on mixed meshes improves total performance

and accuracy of the model in contrast to triangular meshes [27], for quadrilateral cells involving fewer edges than triangular cells. Newly developed methods for postprocessing on unstructured meshes with a rapidly growing community significantly improved the overall performance of such modeling, from preparing the setup through to final plotting. Without the problems related to one- or two-way nesting and with significantly reduced problems at the open boundary (because the meshes are compatible with global models), FESOM-C has become a realistic substitute for existing models when it comes to simulating on regional scales, from particular events to climate simulations.

#### **6. Code Availability**

The version of FESOM-C v.2 used to carry out the simulations reported here can be accessed from https://doi.org/10.5281/zenodo.2085177. The datasets needed for running this setup are available from the corresponding author upon a reasonable request.

**Author Contributions:** I.K., A.A., and V.F. designed experiments. I.K. setup and carried out the experiments. I.K. analyzed and visualized the model as well as the observed data. I.K. wrote the paper with the support of A.A., V.F., and S.D. A.A. analyzed and wrote up the results in the section on "solution convergence on different meshes". A.A. developed the FESOM-C model with the support of I.K., V.F., and S.D. S.H. and N.R. carried out code optimization and parallelization. I.K., A.A., V.F., and S.D. contributed with discussions of the results. K.H.W. helped supervise the project. All authors discussed the results and commented on the paper at all stages. All authors have read and agreed to the published version of the manuscript.

**Funding:** We acknowledge support by the Open Access Publication Funds of Alfred-Wegener-Institut Helmholtz-Zentrum für Polar- und Meeresforschung. This research was partly funded by the state assignment of FASO Russia (theme 0149-2019-0015).

**Acknowledgments:** We would like to thank Semjon Schimanke from SMHI for help in preparing and providing the atmospheric forcing. We would also like to thank the teams behind the COSYNA, ICES, and EMODnet databases for data gathering and maintenance. We thank the administrators of the AWI cluster Ollie, where the final simulations took place, for their continual support and patience. We moreover thank the many people from the Institute of Coastal Research Helmholtz-Zentrum Geesthacht (HZG), where IK began preparing the current manuscript, for their great cooperation and useful discussions, in particular: Holger Brix, for general support at the early stages of preparing the current manuscript; Yoana G. Voynova and Onur Kerimoglu, for discussing the FerryBox data; and Onur Kerimoglu again, for sharing river discharge data. We thank reviewers who have contributed their time and expertise to review manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
