Sensitivity of results of Black Sea level modeling to the choice of boundary conditions on the free surface

This work is devoted to an analysis of the influence of sea surface boundary conditions on results of Black Sea level modeling. Two numerical experiments are carried out by using a z-coordinate model of the Marine Hydrophysical Institute, RAS. A linear kinematic condition on the free surface is used for calculating a level rise in the first experiment and a nonlinear one – in the second experiment. The analysis shows that both experiments reconstruct the Black Sea surface dynamics in accordance with altimetry data. It is found that the sea level near the western boundary of the basin decreases under the nonlinear kinematic condition in cases of storm winds. The greatest differences in the structure of the sea level fields are observed during weakening of the RIM Current in the evolution zones of Sevastopol and Batumi anticyclones. The minimal and maximal values of the level are almost identical in both experiments, but the locations do not match. The discrepancies in the spatial positions of raising/lowering levels are due to a mismatch in the motion phases of anticyclones and mesoscale eddies along the periphery of the RIM Current.


Introduction
Level height variations in coastal zones of the enclosed seas can pose a serious threat to the population and objects of state and private business. As shown in the observational data, catastrophic level rises in the Black Sea and the Azov Sea reached heights of about 7 m [1], and in the Caspian Sea, of more than 20 m [2], which caused flooding of the territory and significant damage. Such extreme situations were induced by meteorological events. On the other hand, over the second part of the twentieth century a trend to increase of the Black Sea level was observed, which was mainly connected with the growth of the surface freshwater flux calculated as the difference between the precipitation and evaporation [3]. The tide-gage and satellite altimetry data of 1992 -2005 also show that the rising rate of the mean sea level is 7.6 ± 0.3 mm/year, which is 2 -3 times higher than the estimates for the previous periods [4]. In addition to a positive trend of the Black Sea level into a basin scale, the spatial heterogeneity of the energy distribution of the level oscillations was also noted. According to [5], the shallow area of the North-Western Shelf exposed stronger variability in the synoptic and mesoscale ranges than the abyssal zones near the Crimean and eastern coasts.
Therefore, accurate prediction of the sea level changes is important from the point of view of risk assessment for the coastal areas. One of the main ways to obtain information about the past, present,  2 and possible future state of ocean systems is numerical modeling. The quality of the diagnosis and prognosis depends on correct consideration of numerous factors. In particular, the accuracy of the sea level reconstruction in the hydrodynamic models is determined on the basis of adequate formulation of boundary conditions on the sea surface. This work is devoted to the analysis of modeling results of the sea level fields for different surface boundary conditions. Section 2 describes the problem formulation. The numerical model being used and the parameters of the experiment are presented in Section 3. A comparison of the modeling results of the Black Sea level fields is made in Section 4. Analysis of the eddy kinetic energy is performed in Section 5. The main conclusions are presented in Section 6.

Problem formulation
There are two main approaches to constructing circulation models. In the first case, the solution of the ocean thermodynamics equations is constructed in terms of a stream function (for example, the MICOM [6]) and in the second one, in terms of a reference ocean level (for example, POM [7], NEMO [8]). The second approach allows one to calibrate models by taking into account in-situ observations of the sea level. Earlier the rigid-lid condition was used in works on modeling circulation in the Black Sea [9]. This condition filters fast inertial-gravity waves. As shown in paper [10], in the shelf areas of the Black Sea the dynamics of the currents during the adaptation period after the initial perturbation differs significantly from the data obtained with the rigid-lid condition. Consequently, inertial-gravity waves make a considerable contribution to the process of adaptation of density fields to the bottom relief and shore orography. Thus, it is necessary to specify more rigid conditions on the sea surface.
The free surface equation is constructed under the assumption that the kinematic condition is satisfied, indicating that the free surface is impenetrable for the liquid: Since the vertical velocity of the surface motion is defined as the total time derivative from z, the free surface equation has the following form for case of absence of a mass flux from the atmosphere: where u, v, w are the current velocity components on the free surface; ζ is the sea level elevation. If the ratio ζ/H << 1 (where H is depth), the third and fourth terms in the left-hand side of equation (1) are neglected (the linear condition). The maximum depth of the Black Sea exceeds 2 km, while more than 20% of the sea area is occupied by domains with depths less than 100 m. The North-Western part of the sea is an extensive shelf zone with average depths of about 20 m. As shown in book [1], extreme level rises can reach here several meters and, therefore, the condition ζ/H << 1 becomes incorrect in some cases. The aim of this work is to estimate the accuracy of the modeling results of the Black Sea level fields by using the linear and nonlinear kinematic conditions on the free surface.

Model and experimental setting
An eddy-resolving z-model of the Black Sea circulation developed at the Marine Hydrophysical Institute of the Russian Academy of Sciences (MHI RAS) [11] was used in our study. The model is based on the primitive equations of ocean taking into account the approximations of Boussinesq, hydrostatics, and incompressibility of the ocean water: Here f is the Coriolis parameter; g is the gravity acceleration; p is the pressure; T is the temperature; S is the salinity; ρ is the density; ρ 0 = 1 g cm -3 ; ρ 1 is the mean density in the first model horizon; ν V and ν H are the vertical and horizontal viscosity coefficients; κ V and κ Н are the vertical and horizontal diffusion coefficients. The density nonlinearly depends on the temperature and salinity (equation (7)). The density change due to pressure was neglected according to [12]. The wind stress (τ x , τ y ), the heat flux Q, the precipitation Pr , and the evaporation Ev are set on the sea surface as boundary conditions: , , where S 0 is the surface salinity; cl S is the surface climatic salinity; and  is the relaxation parameter. The relaxation parameter  is equal to 2/1728 m day -1 and it denotes the salinization velocity. We set the non-flow condition on the solid boundary. The water exchange through the Bosporus and the Kerch Strait and the runoff of the main Black Sea rivers are taken into account by the Dirichlet boundary condition. The impermeability condition and the absence of heat and salt fluxes are set on the bottom, and the bottom friction is not considered. The slip condition is satisfied, and the normal derivatives for temperature and salinity are zero on the solid lateral boundaries.
The model is implemented on a C-grid [13]. The advective terms in equations (5), (6) are approximated by the total variation diminishing (TVD) schemes [14]. Biharmonic operators for the calculation of the horizontal viscosity and diffusion in equations (2), (3), (5), (6) were used to provide filtering of the computational noise and stabilizing of the numerical solution. The vertical turbulent mixing is parameterized by the Mellor-Yamada turbulent closure model of level 2.5 [15]. Full description of the model, mathematical assumptions, and numerical approximations are presented in [11].
We carried out two numerical experiments on Black Sea circulation modeling taking into account real atmospheric forcing. The difference between the experiments was in a form of the free surface equation. The linear kinematic condition (equation (8)) was set in the first experiment, and we used the nonlinear kinematic condition (equation (9)) in the second one: The experiments were performed on a uniform grid with a horizontal step of 1.6 km, 27 z-horizons were set vertically, and the time step was 96 sec. The biharmonic operator coefficients in terms describing the horizontal turbulence in equations (2), (3), (5), (6) are ν H = 10 16 cm 4 s -1 and κ H = 10 16 cm 4 s -1 . The bathymetry was constructed by using the Black Sea depths from the MHI oceanographic data bank. The sea level, temperature, salinity, and horizontal velocity fields were set at the initial moment of time. These data correspond to January 1 of a climatic year. The data for 2011 obtained by the SKIRON/Eta model with a horizontal resolution of 0.1° in latitude and longitude were used as atmospheric forcing. Observation data are not assimilated in the model. All boundary and initial fields were linearly interpolated to the model grid nodes. Both experiments were carried out for one year. The daily-mean fields of all hydrophysical parameters were returned as output data. Validation of the modeling results by in-situ data is presented in [16].

Modeling results
Let us consider the calculated sea level elevation ζ for two experiments. Subscripts L and N in the text below will be used to denote the results of experiment 1 (linear case) and experiment 2 (nonlinear case), respectively. The difference Δζ = ζ L -ζ N was calculated at each grid node. Maps of the spatial distribution Δζ were created to quantify the estimation of the discrepancy between the results of the two calculations. The analysis of maps of the level difference fields for each day showed that the values Δζ did not exceed an average of 0.6 cm during the first three months after the beginning of simulation. A fairly smooth structure of the currents field was observed in the upper layer in the winter of 2011. Therefore, the sea level fields were reconstructed almost identically in both experiments. The accuracy of the modeling results was assessed based on data of the operational system for diagnosis and forecast of hydrophysical characteristics of the Black Sea [17]. Note that these materials were calculated with the assimilation of available observational data (including altimetry). To estimate the correspondence of our results to reanalysis data, the differences between the year-mean sea level obtained in [17] and the year-mean values of ζ L and ζ N were compared. Before that the experimental outputs were interpolated to the reanalysis model grid (resolution: 5 km). Spatial distributions of the difference are presented in Figure 1. The largest deviations (3.4 and 2.7 cm for experiments 1 and 2, respectively) were observed in the abyssal part of the sea, and they are associated with the initial conditions in the experiments. Recall that we used, as the starting point, climatic hydrophysical fields which are smoother and lower than real data. As we see, the maximum absolute magnitude of the sea level difference decreased for experiment 2, especially in the northern and south-western parts. A good agreement between the daily-mean sea level fields and reanalysis data was obtained for situations when the level height increased two and more times compared with the year-mean maximum absolute values.    Significant differences in the level field structure were revealed approximately 7 months after the start of the simulation. Spatial distributions of the time-averaged Δζ calculated for January-July 2011 and August-December 2011 are presented in Figure 3. One can see that extremum values of the timemean Δζ increased by an order in the second part of the year compared with the first half. For the daily-mean magnitudes the maximum Δζ did not exceed 10% of the maximum values ζ L and ζ N in the first half of the year. However, occasionally in the second part for both experiments Δζ reached 15 cm and was commensurate with the maximum level height.
The Black Sea circulation is characterized by a sub-basin cyclonic gyre, the so-called RIM Current propagating along the continental slope. The RIM Current weakens in the warm season and intensifies in autumn-winter. Anticyclonic eddies are mainly formed between the RIM Current and the shore. The most intensive and long-lived eddies in the Black Sea are the Sevastopol and Batumi anticyclonic eddies. Their size is of the order of 100 km and their life-time varies from 4 to 8 months. Relatively the Rossby deformation radius for the Black Sea (about 25 km), these eddies are mesoscale and quasigeostrophic (the Rossby number << 1) in accordance with a classification of the oceanic processes of [18]. Figure 3 (b) illustrates that the greatest differences in the level are located in the center of the sea eastern part and in domains of the Sevastopol and Batumi anticyclones. A comparison of the daily-mean sea level and surface current fields revealed that the displacement phases of the mesoscale eddies moving along the periphery of the RIM Current did not coincide with the experimental results. It was found that the Sevastopol anticyclonic eddy began to move to the south-west about 40 days earlier in experiment 2 than by data of experiment 1 in the second part of the year. A cyclonic eddy located to the north-east of the Batumi anticyclone was more intense and occupied a large area according to the results of the first experiment.  In experiment 2, the Sevastopol anticyclone did not disappear completely after the storm, but only slightly weakened on October 27 -30, in contrast to experiment 1. It shifted along the core of the RIM Current to the south-west and dissipated in mid-November. A new anticyclone generated to the southwest of Sevastopol in October occupied its place in early November. The eddy increased in size and moved along the RIM Current till the end of the year. During December a chain of anticyclonic eddies (rather than individual cyclones, as in experiment 1) was formed near the Crimean coast. This structure was moving to the west not enhancing the Sevastopol anticyclone. It was found that the Sevastopol anticyclone was less intense than in experiment 1 and its size varied in the range of 100 -150 km.
The Batumi anticyclone had almost the same size and intensity by the beginning of the storm on October 17 -18 in experiments 1 and 2. It dissipated a few days after the storm, the RIM Current being pressed to the shore in both cases. Since November, the formation of coastal anticyclones in the south-eastern part was reconstructed in the results of the numerical experiments differently. The coastal eddies were more intense and larger in area for experiment 2. During December, the number of anticyclonic formations between the RIM Current and the shore in both experiments was different. Their number remained unchanged and was equal to two in the second experiment, whereas the number of eddies varied from two to four and their sizes were smaller in the first experiment. Cyclonic eddies were formed between the anticyclonic ones in December, but they were more extensive in experiment 2.

Energetics
As shown above, the dynamics of mesoscale eddies changed when applying the nonlinear kinematic condition on the free surface. The mesoscale variability of circulation is usually associated with the eddy kinetic energy determined by the fluctuations of the velocity of currents from the mean value [19]. The Lorenz Energy Cycle analysis [20] is generally used to study the mechanisms of circulation variability. The energy cycle is defined by the balance of four components: the mean kinetic energy (MKE), the eddy kinetic energy (EKE), the available mean potential energy (MPE), and the available eddy potential energy (EPE). To investigate the possible reasons of distinctions in the eddy displacement phase, the EKE variability for both experiments was considered. The potential energy components were not investigated here. The year-mean (time-averaged over 2011) and time-varying (deviations from the year-mean for each day of the year) zonal and meridional velocity components were calculated. Based on these data, the volume-integrated values of the EKE and the rate of energy conversion from the EKE to the MKE through barotropic instability (denoted by EKE and С(EКЕ, MКE), respectively) were calculated for each day by using the following formulas: where  V dV is the integral over the basin volume; the apostrophe denotes the deviation from the yearmean value. If С(EКЕ, MКE) is positive, the energy is converted from the EKE to the MKE. The yearmean value of the volume-integrated kinetic energy (МКЕ) was equal to 0.563 PJ (petajoules = 10 15 joules) for both experiments. Also, it was found that the mean value of EKE L was equal to 0.228 PJ and the mean value of EKE N was equal to 0.217 PJ. Thus, in 2011 the year-mean contribution of the EKE to the kinetic energy was 40.5% for experiment 1 and 38.5% for experiment 2.
Changes in the wind stress maximum over the Black Sea, EKE and С(EКЕ, MКE) during 2011 for both experiments are shown in Figure 4. One can see in Figure 4

Conclusions
A study of results of Black Sea level modeling calculated by using linear and nonlinear kinematic conditions on the free surface was performed. Experimental data were obtained by an eddy-resolving numerical MHI-model with a resolution of 1.6 km with real atmospheric forcing for 2011. The dailymean and year-mean values of the sea level were analyzed. The differences between the year-mean reanalysis sea level and year-mean values of ζ L and ζ N were compared to estimate the correspondence of the modeling results with the reanalysis data (with altimetry assimilation). It has been revealed that the deviation of the sea level height from the reanalysis data decreased in the experiment with the nonlinear kinematic condition. It has been found that the daily-mean height of the sea level near the western boundary of the basin was decreased by using the nonlinear kinematic condition for cases of significant level rises (more than 25 cm). The reduction of the maximum level height in the coastal zones can be explained by the fact that the account of all terms in the free surface equation allows us to describe a wider class of level oscillations. As a result, in the second experiment the eddy kinetic energy was decreased due to a more correct redistribution between the movements on different spatial and temporal scales. To construe the types of level oscillations reconstructed by using the nonlinear kinematic condition, we plan to perform a spectral analysis of the modeling results.
The comparison of the experimental results demonstrated that the greatest differences in the structure of the sea level fields were observed during the period of weakening of the RIM Current (the warm season) in the evolution areas of Sevastopol and Batumi anticyclones. Note that the minimum and maximum values of level heights are almost identical in both experiments, but their locations do not match each other. The discrepancies in the spatial positions of raising/lowering levels are due to a mismatch in the displacement phases of anticyclonic eddies along the periphery of the RIM Current. The Sevastopol and Batumi anticyclones are components of mesoscale variability in the Black Sea. In the second part of 2011, changes in the dynamics of mesoscale eddies were observed. This process was displayed in the behavior of the volume-integrated eddy kinetic energy when the discrepancy between EKE L and EKE N became significant after a storm in July and reached a maximum after a storm in October 2011. The restructuring of the Black Sea mesoscale circulation in 2011 in the second experiment is associated with the decrease in the eddy kinetic energy due to a reduction in the inflow from the mean current kinetic energy. It is probably related to a redistribution of kinetic energy towards the submesoscale processes.
The above-obtained results show that the account of nonlinear terms in the equation for calculating the free surface height led to a qualitative change in the response of the Black Sea current field to an intense wind impact in 2011. Thus, the deviations of the level heights from the reanalysis data were reduced. These conclusions should be taken into account for improving the accuracy of forecasting the hydrophysical fields in the Black Sea, especially in the submesoscale ranges of variability. Now the operational system for diagnosis and forecast uses a version of the MHI-model where the level height