Introduction

Major Baltic inflows (MBIs) represent the irregular introduction of extremely large volumes of the North Sea waters, from 90 to 258 km³, into the Baltic Sea, lasting 6–29 days, which penetrate into the deep-water areas of the Baltic Proper (Fig. 1), exerting a beneficial effect on the ecological state of this sea [1–7]. Weak inflows of the North Sea waters of 10–20 km³ in volume appear constantly but the penetration of these waters into the Baltic is most often limited to the Arkona Basin (Fig. 1). MBIs are a relatively rare phenomenon observed until the early 1980s from 1–2 times a year to once every few years [4]. Spreading far into the open part of the Baltic Sea, highly saline and oxygen-rich waters of large inflows renew the bottom and deep-water Baltic masses exposed to hypoxic conditions  [2, 4, 8]. Observations show that after 1983, the number of MBIs decreased significantly and the interval (stagnation period [8]) between them increased greatly and amounted to 10–11 years [4, 6, 7–9]. Physical mechanisms for the increase in stagnation periods remain unexplored to date. The last major inflow occurred in December 2014 [5], after which no new MBIs have been described in the scientific literature.

MBIs can be considered as an extreme water exchange component between the North and Baltic seas. For example, according to K. Wyrtki [10] and H. Fischer and W. Matthäus [3], about 200–225 km³ of the Kattegat waters passed through the Danish straits during the MBI in November – December 1951, which amounted to approximately 40% of the annual norm.

The accumulated information on the variability of hydrometeorological processes during MBI allowed scientists to identify four periods in the process of its formation: outflow period, precursory period, main inflow period and post inflow period [4, 5, 8].

The outflow period starts when eastern winds blow over Northwestern Europe, which contributes to the water outflow from the Baltic into the North Sea and its level decrease. This period is very important for the future formation of MBIs since the longer and more intense the water outflow from the Baltic is, the more its level will decrease and the greater the level gradient between the Kattegat and the Southwestern Baltic will form before the beginning of MBIs. The intensity of MBIs depends largely on this gradient [4, 5, 8].

F i g.  1. Bathymetry of the North and Baltic seas (black square indicates the area of the Southwestern Baltic and the Kattegat) (a), enlarged image of the selected area (b). Designations: asterisks show location of the Darss Sill (DS) and Arkona (A) automatic stations; AB is the Arkona Basin, BD is the Bornholm Deep, GD is the Gotland Deep, LD is the Landsort Deep, LB is the Little Belt Strait, GB is the Great Belt Strait; S is the Sound Strait, SK is the Skagerrak, MB is Mecklenburg Bay, KB is Kiel Bay, GG is the Gulf of Gdansk

In the precursory period, the synoptic situation changes: the east wind weakens and begins to change its direction to the west one. This leads to the sea level elevation in the Kattegat, gradually approaching the level in the Southwestern Baltic [4, 5, 9].

The main inflow period occurs when the North Sea level elevation, which began in the previous period, reaches a critical value, at which the level gradient becomes directed from the Kattegat to the Southwestern Baltic and continues to grow under the influence of strong west winds, the duration of which reaches 2–3 weeks. At this time, the difference in level between the Kattegat Strait and the Southwestern Baltic (Fig. 1, b) can reach 1.0–1.7 m [11]. As a result, large masses of the highly saline and oxygen-rich Kattegat waters enter the Baltic Sea, which, in turn, leads to a further decrease in the North Sea level and an increase in the Baltic Sea [4, 5, 8].

The post inflow period begins when the west winds weaken and the North Sea waters cease to accumulate in the Danish straits. Since the Baltic Sea level is elevated relative to the North Sea level, the water outflow from the Baltic Sea begins and its level drops to a level close to its average value [4, 7, 8].

Mathematical modeling of water exchange and oceanographic conditions in the North and Baltic Seas is a complex task for two main reasons. Firstly, the oceanographic regimes of these seas are very different. The North Sea is a shallow (except for the Norwegian Trench) (Fig. 1, a), weakly stratified marine basin with intense tidal dynamics and vertical mixing, relatively freely communicating with the ocean, so its salinity is close to oceanic. On the contrary, the Baltic Sea is an almost completely closed brackish marine basin with very weak tidal dynamics and sharp stratification, limiting vertical mixing between surface and deep-water masses. Secondly, due to the narrowness and shallowness of the Danish straits connecting the North and Baltic Seas (the Sound, the Great Belt and the Little Belt) (Fig. 1, b), which have a complex morphometry of the coastline and bottom topography. The minimum width of the Sound is less than 5 km and its smallest depth is 8 m; for the Great Belt Strait, these parameters are 3.7 km and more than 20 m; for the Little Belt Strait, they are 0.8 km and 12 m, respectively [8, 12, 13].

Such characteristics of the Danish straits require the use of a grid domain with cells of significantly smaller dimensions than the smallest width of the straits in numerical modeling for correct water flow simulation in these straits as well as the features of stratification and current structure. Processing power of modern computers does not make it possible to use uniform grids with such a high spatial resolution for modeling not only the combined water area of the North and Baltic seas, but also the Baltic Sea alone. To solve this problem, scientists expanded the Danish straits artificially when modeling the oceanographic conditions of the Baltic Sea, adjusting their width to the spatial resolution of the grid domain used in the model [14–16]. Such a procedure, with an unchanged depth, led to a change in the cross-sectional area of the straits. Therefore, the depth of the straits was reduced to maintain the cross-sectional area. Both changes in the morphometry of the straits lead to changes in stratification, current structure and salt transport volume in the Danish straits [12].

An important length scale that should be taken into account when modeling oceanographic fields for correct resolution of mesoscale eddies, upwellings [17] and structure of narrow jets caused by the dynamics of gravity currents in the Southwestern Baltic [18, 19] is the baroclinic Rossby radius of deformation. According to estimates by various researchers, its largest values (7–9 km) were observed in the Bornholm Basin and deep-water areas of the Baltic Proper, with the smallest (1–2 km) ones in the sea shallow areas with depths less than 50 m [20–23]. In this regard, models with nested grids have been used to improve the spatial resolution in numerical experiments. For example, in [12], the model domain with a nested uniform grid had a spatial resolution of 900 m and included the waters of the Kattegat, the Danish straits and the Arkona and Bornholm basins of the Baltic (Fig. 1, b). One of the liquid boundaries of the model was located in the north of the Kattegat and the other one – in the east of the Bornholm Basin [12]. However, such models do not allow studying the propagation routes of the MBI waters in other Baltic Sea areas.

In our opinion, models with unstructured grids with the highest condensation (detailing) in the area of the Danish straits are more promising for the MBI studying, which allows for a more accurate description of the structure of currents, stratification of water masses and salt transport through narrow and shallow straits. In [24], a joint model of the North and Baltic Seas with a mixed triangular-quadrangular unstructured grid was used, which made it possible to achieve a nominal spatial resolution of 200 m in the Danish straits. Comparison of the modeling results with data from tide gauge measurements of sea level and measurements of temperature and salinity at the Fehmarn Belt and Arkona stationary stations showed good agreement overall, although in some areas of the compared series, the discrepancies between the measured and calculated values reached 30–50 cm for sea level, 3–5°C for temperature and 2–3‰ for salinity [24].

The main purpose of the paper is to estimate the capabilities of numerical hydrodynamic modeling of MBIs using a 3D baroclinic model of the North and Baltic Seas, which has a spherical grid area with detailing in the Danish straits, and, based on the modeling results, to study the structure and propagation routes of the transformed North Sea water flows in the Baltic Sea after the MBI in December 2014.

Data and methods

Model description

The Institute Numerical Mathematics Ocean Model (INMOM) of ocean and sea circulation developed at Marchuk Institute of Numerical Mathematics of RAS [25, 26] was chosen as the basic model to describe oceanological processes in the Baltic and North Sea systems during the 2014 MBI.

The INMOM is based on a complete system of nonlinear primitive equations of ocean hydrodynamics in spherical coordinates in the hydrostatic and Boussinesq approximations. Dimensionless quantity σ=(z-ζ)/(H-ζ)  is used as a vertical coordinate, where z is usual vertical coordinate; ζ=ζλ,φ,t  is deviation of the sea level from the undisturbed surface as a function of longitude λ, latitude φ and time t; H=Hλ,φ is sea depth. The number of vertical sigma layers in the model is 20.

The predictive variables of the model are horizontal components of the velocity vector, potential temperature, salinity and ocean level deviation from the undisturbed surface. To calculate the density, a specially designed for numerical models [27] equation of state that takes into account the compressibility of seawater is used.

The coefficients of vertical turbulent diffusion and viscosity were selected according to the Pacanowski – Philander parameterization [28]. The coefficient of turbulent diffusion varied from 1 to 50 cm²/s and the coefficient of turbulent viscosity varied from 1 to 250 cm²/s. Horizontal turbulent diffusion and viscosity were described by the usual Laplacian with coefficients ν  = (3–8)×10⁴ cm²/s. Bottom friction was specified by a quadratic equation with coefficient C_D = 2.5×10⁻⁴.

The model includes a sea ice thermodynamics block [29] consisting of three modules. The thermodynamics module describes freezing, ice melting and snowfall. The ice dynamics module calculates its drift velocities. The ice transport module is used to calculate the ice and snow cover evolution due to the drift [30].

The model uses a spherical grid with two poles, one of which is located on the Jutland Peninsula (Denmark) and the other – in the very south of Sweden (Fig. 2). The spatial resolution of the grid area nodes in the area of the Danish straits is about 300–700 m and increases proportionally to 6–12 km with distance from the straits towards the outskirts of two seas.

F i g.  2. Grid area of the model. Red dots indicate model liquid boundaries, black circles – grid area poles

For this model version, bathymetry from GEBCO 2015  was combined. When preparing the model bathymetry, depth values were interpolated into grid nodes and smoothed with a Gaussian filter to eliminate their sharp differences, which increases significantly the stability of calculations during modeling.

Initial and boundary conditions

The initial conditions were monthly average water temperature and salinity data for January 2014 with a vertical resolution of 5 m and a spatial resolution of 4.5 × 9 km from the GLORYS12V1 ocean reanalysis  (available at: http://marine.copernicus.eu).

For the boundary conditions on the sea surface in the INMOM atmospheric module, such meteorological characteristics as air temperature and humidity at a height of 2 m, pressure at sea level, wind speed at a level of 10 m, incident shortwave and longwave radiation, atmospheric precipitation were specified with a discreteness of 3 hours, a spatial step of 0.25° and duration from January 2014 to December 2015 obtained from the ERA5 reanalysis .

At the liquid boundaries of the North Sea (Fig. 2), the average monthly values of water temperature and salinity observed from January 2014 to December 2015 as well as the amplitudes and phases of the oscillations of the level and currents of eight main tidal harmonics (M2, S2, N2, K2, K1, O1, P1, M4) taken from TPXO9 global tidal model (available at: https://www.tpxo.net/global) were specified.

On the solid sections of the lateral boundary, the heat and salt fluxes were set equal to zero and the no flow and free sliding conditions were used for the current velocity.

Model calculations were carried out from 1 January 2014 to 31 December 2015, with the average results being derived for each hour.

Verification of the model and comparison of the modeling results with the data of regional reanalysis of hydrophysical fields

To verify the model, we used the data from contact measurements of salinity and currents at different horizons of stationary automatic stations Darss Sill and Arkona located in the Southwestern Baltic at 21 and 45 m depths, respectively (Fig. 1, b). Observations of salinity at the Darss Sill station are made at 2, 5, 7, 12, 17 and 19 m horizons and at the Arkona station – 2, 5, 7, 16, 25, 33, 40 and 43 m. Current velocity and direction are measured with Doppler acoustic profilers at these stations.

The INMOM modeling results were compared with instrumental measurements as well as with salinity and current change data obtained using the regional reanalysis of the Baltic Sea Physics Analysis and Forecast (BSPAF) hydrophysical fields based on the numerical implementation of the Nucleus for European Modeling of the Ocean (NEMO) 3.6 hydrodynamic model [31, 32] for the Baltic Sea conditions. This model uses a procedure for contact and satellite information assimilation based on the algorithm of one of the varieties of the Kalman filter (local singular evolutive interpolated Kalman (LSEIK) filter) [33]. Satellite data on surface water temperature provided by the Swedish Meteorological and Hydrological Institute (SMHI) ice service as well as in situ T and S measurements from the ICES database (available at: http://www.ices.dk) were used as assimilated variables in the NEMO 3.6 model. The NEMO 3.6 model used meteorological data computed with the ECMWF ERA5 atmospheric model to set the sea surface boundary conditions. The BSPAF regional marine reanalysis data are daily averaged, they have a horizontal resolution of 3.9 km and 56 vertical horizons (layer thickness varies with depth from 3 to 22 m) and cover the 1993–2022 period.

To compare the salinity changes measured and calculated with the INMOM model and the BSPAF reanalysis data at different depths, mathematical expectations of salinity series m_s, their standard deviations (SD) s_s, minimum S_min and maximum S_max values and correlation coefficient R_ss between the measured and model salinity values were estimated. The accuracy of the salinity values calculated using the INMOM and NEMO 3.6 models (BSPAF reanalysis) was estimated using accuracy criterion P_a, which shows the number of salinity values calculated using the models that fall within range < 0,674s, where s is SD of the salinity values measured at the Darss Sill and Arkona stations.

To compare the measured and model values of currents, the following statistical characteristics of the variability of the velocity and direction of currents were estimated using the vector-algebraic method of analysis of random processes.

1) mathematical expectation of vector process m_v (module direction m_V  and direction a_m);

2) linear invariant of the SD tensor [I₁(0)]^0,5, where I₁(0) = l₁(0) + l₂(0) is linear invariant of the vector process dispersion tensor determined through the half-lengths of principal axes l₁(0) и l₂(0) of the dispersion ellipse and orientation a° of its major axis relative to the geographic coordinate system:

λ1,2(0)=12Dvv+Duu±Dvv-Duu2+Dvu+Duv2 ,

where Dvv,Duu  are dispersions of the vector process components;

3) stability of currents  ½m_v ½, where ½m_v ½ is module of the mathematical expectation of a vector process. When r > 1, the intensity of oscillatory motions in the flow prevails over the intensity of the average transfer, i.e. the current is unstable; when r < 1, currents become more stable;

4) two invariants of the normalized cross-correlation tensor function between the currents measured at the Darss Sill station and calculated using the INMOM model and BSPAF data: linear invariant I₁^VU(t) and rotation indicator D^VU(t). Linear invariant I₁^VU(t) is equal to the trace of the matrix of correlation tensor function K_VU(t), two vector processes V(t) and U(t) and characterizes the integral of the intensities of collinear changes in vector processes V(t) and U(t):

K_VU(t) = &Kv1u1(τ),Kv1u2(τ)&Kv2u1(τ),Kv2u2(τ) ,

where t is time shift; v₁ is component of vector process V(t) on the parallel; v₂ is component of vector process V(t) on the meridian; u₁ is component of vector process U(t) on the parallel; u₂ is the component of vector process U(t) on the meridian.

Rotation indicator D^VU(t) is equal to the difference of the non-diagonal components of the matrix of correlation tensor function K_VU(t) and characterizes the integral of orthogonal changes in processes V(t) and U(t); when D^VU(t) > 0, process U(t) is rotated on average relative to process V(t) over a given time period clockwise and counterclockwise when D^VU(t) < 0.

Then the total correlation coefficient was calculated:

R_VU(t) = [I1VU(t)]2+[DVU(t)]2.

In addition, the maximum modules of current velocity V _max were estimated.

Flow rates of currents Q through the Danish straits during the 2014 MBI formation were estimated based on the current velocity vectors (V) calculated by the INMOM model at different horizons along three sections crossing the straits (see Fig. 1, b), using the following formula:

Q=i=1nz=1mVS,                                                       (1)

where n is number of i-cells on the section; m is number of z-horizons in the given cell; V  is meridional component of the current velocity in a grid cell at horizon z; S is cross-sectional area of the cell determined as product of layer thickness (∆z) and distance between adjacent nodes of the grid region of model (∆i ), i.e. S=∆z×∆i .

To study the propagation routes of transformed North Sea waters after the MBI, two methods were used. The first method made it possible to construct two oceanographic sections along a system of interconnected deep-water basins of the marine relief. Their location was determined based on published information on the migration routes of the salty North Sea waters during the MBI in the Baltic Sea [4, 5]. Using the modeling data, diagrams of the temporal variability of salinity in the bottom layer were constructed on these sections. In the second case, the Lagrangian method was used.  Its detailed description is given in [34]. Within the framework of this method, 5000 passive markers were placed daily from 1 November to 31 December 2014 on a segment along the boundary passing north of the Danish straits along a line with coordinates 56.6°N, 10.85°E – 56.6°N, 11°E (Fig. 1, b). Based on the current velocity vector fields calculated using the INMOM model, each marker trajectory was calculated for a period of one year (until 31 December 2015).

Lagrangian trajectories were calculated using advection equation

dλdt=uλ, φ, t,

dφdt=vλ, φ, t,

where u and v are angular components of the current velocity calculated using the INMOM model in the penultimate s-layer in depth; φ and λ denote latitude and longitude, respectively. Angular velocities are used to simplify the equation of motion on a sphere. Velocity values inside the grid cells were calculated using bicubic interpolation in space and third-degree Lagrange polynomial interpolation in time. When modeling the Lagrangian trajectories, the coordinates of the passive markers were recorded with the 1 h time resolution.

Results and discussion

Comparison of salinity and current values measured and calculated with the INMOM model and BSPAF reanalysis

Figs. 3 and 4 show changes in salinity values obtained during measurements at different horizons of the Darss Sill and Arkona automatic stations (see Fig. 1, b for the location of stations) based on the INMOM modeling results and the BSPAF regional reanalysis data for 1 November – 31 December 2014. Table 1 shows the statistical estimates of the measured and model salinity values. It is evident that the December 2014 MBI is modeled both by the BSPAF regional reanalysis data and by the INMOM modeling results as an anomalously large increase in salinity from the bottom horizons to the sea surface (Figs. 3 and 4). At the same time, the BSPAF reanalysis data, unlike the INMOM model, did not reproduce two weak inflows of the Kattegat waters that occurred on 22 and 26 November 2014 (Fig. 3). The correlation coefficients (R_ss) between the measured and model (INMOM and BSPAF) salinity series at different horizons are high and vary from 0.70 to 0.98 near the Darss Sill station and from 0.67 to 0.98 near the Arkona station (Table 1). This result indicates that the reanalysis data and the INMOM model describe adequately the main features of salinity changes during the MBI in the Southwestern Baltic, although the values of correlation coefficient R_ss between the measurement results and the INMOM data for three upper horizons near the Darss Sill station are noticeably higher than those of BSPAF, while they are close for three lower horizons. On the contrary, for the Arkona station area, the R_ss values at three upper horizons are lower for INMOM compared to BSPAF and the R_ss values at three lower horizons for INMOM are higher than those of BSPAF (Table 1).

F i g.  3. Water salinity at the Darss Sill station based on the measurement (a), BSPAF reanalysis (b) and INMOM modeling (c) data from 1 November to 31 December 2014

F i g.  4. Water salinity at the Arkona station based on measurement (a), BSPAF reanalysis (b) and INMOM modeling (c) data from 1 November to 31 December 2014

The values of the mathematical expectation of salinity changes during the BSPAF formation and propagation period estimated by INMOM in the Darss Sill station area are overestimated by 9–21% relative to the measured values at almost all horizons and by 4–22% in the area of the Arkona station. The exception is the 40 m horizon where the INMOM results in the Arkona area showed an underestimation of the mathematical expectation of salinity changes by 7% (Table 1). In contrast to the INMOM results, discrepancies between measured and calculated values of the mathematical expectation of salinity changes based on the BSPAF reanalysis data are generally significantly smaller and vary from 0.3 to 7% in the Darss Sill station area and from 2 to 17% near the Arkona station (Table 1). Only at the 25 and 33 m horizons in the Arkona area, the excess of the mathematical expectation values according to the BSPAF data relative to the measurements is noticeably greater than according to the INMOM data (Table 1).

T a b l e  1

Statistical estimates of daily average seawater salinity at different horizons based on the measurements at the Darss Sill (DS) and Arkona (A) stations as well as the INMOM modeling and BSPAF reanalysis data for the period from 1 November to 31 December 2014

N o t e: m_s is average value; s_s is SD; S_min, S_max are minimum and maximum salinity values; R_ss is correlation coefficient between the measured and modeled salinity values; P_a is accuracy criterion for the salinity values calculated by the models

Discrepancies between the SD values of salinity according to the INMOM model and results of measurements in the Darss Sill station area in the upper 2–7 m layer are small and do not exceed ±5…7% (Table 1). However, deeper down, the INMOM results show SD values that are underestimated by 19–23%. On the contrary, SD estimates according to the BSPAF reanalysis data show values that are overestimated by 19–33% at all horizons (Table 1).

The estimates of salinity SD measured at the Arkona station at the 2–7 m upper horizons are very small (0.51–0.54‰), which is 6.5–7.5 times less than the estimates of salinity SD based on the measurements at the Darss Sill station (Table 1). At these horizons, SD of salinity obtained from the modeling and reanalysis results shows overestimated values: 0.85–0.96‰ for INMOM and 1.04–1.10‰ for BSPAF. At a depth over 7 m, the estimates of salinity SD based on the measurements at the Arkona station increase significantly (by 4–8 times). Here, in the 16–33 m layer, the estimates of salinity SD obtained from the INMOM modeling results are underestimated by 10–28% relative to the measurement data and only at the bottom horizon of 44 m, they are overestimated by 10% (Table 1). The estimates of salinity SD obtained from the BSPAF reanalysis results are overestimated everywhere at depths from 16 to 40 m. They are overestimated most of all at the 16 m (38%) and 25 m (23%) horizons and not significantly at the 33 m (2%) and 40 m (5%) horizons (Table 1).

Comparison of the measured and model-calculated salinity minimum values shows that according to the INMOM results, they are always greater than their measured values at the Darss Sill and Arkona stations in all cases. Moreover, these discrepancies with the measured values increase from the surface where they do not exceed 1–8% to the bottom horizons where they reach 19–45%.

The discrepancies between the measured and BSPAF reanalysis-calculated salinity minimum values in the areas of the Darss Sill and Arkona stations are noticeably smaller and do not exceed ±9…17% (Table 1).

Comparison of the model-calculated estimates of salinity maxima (S_max) with their measured values at the Darss Sill and Arkona stations during the MBI formation and spread period shows that they exceed the measured values almost in all cases (Table 1). In the Darss Sill station area, the S_max model values based on the INMOM results are 1–9% higher than the measured values while they are significantly higher according to the BSPAF reanalysis data, amounting to 7–16% (Table 1).

For the Arkona station area, the S_max model estimates for INMOM exceed the measured values by 18–27% while those obtained from BSPAF reanalysis data are greater than the measured values by 14–19% (Table 1). At greater depths (16–40 m), discrepancies between the S_max estimates calculated from INMOM and its measured values are noticeably smaller and vary from −4 to +8%. The S_max estimates obtained from BSPAF reanalysis data exceed its measured values by 6–24% (Table 1).

Estimates of the P_a accuracy criterion show that the INMOM model simulates salinity changes in the Southwestern Baltic in general better than the BSPAF reanalysis (Table 1). In the Darss Sill area, 43 to 57% of the INMOM salinity estimates fall within the range of measured values less than 0.674s while only 16 to 25% of the values from the BSPAF reanalysis fall within this range (Table 1). For the Arkona area, the estimates of the P_a accuracy criterion from the INMOM modeling results vary from 33 to 82% while from the BSPAF reanalysis they do not exceed 7–49% (Table 1).

T a b l e  2

Statistical estimates of current velocity variability at different horizons (H) at the Darss Sill station (DS) based on the measurement, BSPAF reanalysis and INMOM modeling data for the period from 1 November to 31 December 2014