Numerical Modeling of Marine Circulation with 4D Variational Data Assimilation

: The technology is presented for modeling and prediction of marine hydrophysical ﬁelds based on the 4D variational data assimilation technique developed at the Marchuk Institute of Numerical Mathematics, Russian Academy of Sciences (INM RAS). The technology is based on solving equations of marine hydrodynamics using multicomponent splitting, thereby solving an optimality system that includes adjoint equations and covariance matrices of observation errors. The hydrodynamic model is described by primitive equations in the sigma-coordinate system, which is solved by ﬁnite-difference methods. The technology includes original algorithms for solving the problems of variational data assimilation using modern iterative processes with a special choice of iterative parameters. The methods and technology are illustrated by the example of solving the problem of circulation of the Baltic Sea with 4D variational data assimilation of sea surface temperature information.


Introduction
Numerical models of the dynamics of oceans and seas are important components of modern technologies for monitoring and forecasting the global hydrosphere, oceans and inland and marginal seas. They allow one to simulate complex hydrodynamic processes, and to study the fine structure of hydrophysical fields and their spatio-temporal variability [1][2][3][4].
In modern studies of oceanic and marine circulations, three key factors are distinguished for improving the model performance and accuracy of numerical forecasting. The first is to increase the spatial resolution of numerical models. It allows one to explicitly describe the mesoscale eddies which are often highly energetic, to improve and detail the qualitative structure of the simulated hydrodynamic fields [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19]. The second, one of the most interesting points, towards which scientific modeling research is moving, is that of coupled models with a strong multidisciplinary component (atmosphere, ocean, waves, aerosol, biological component, etc.). Coupling plays an important role both for very small space-time scales (order of km and hours) and for very large scales (synoptic scale for short term and climatological approaches) [20][21][22][23]. These questions are not discussed in our article.
The third factor is related to the development of observational systems and numerical techniques for data assimilation. This allows one to link the model trajectory to real data and increase the accuracy of the forecast. Especially in coastal areas, with a greater presence of data in situ, assimilation techniques can play a fundamental role in overcoming the shortcomings of satellite data.
Currently, there is an increasing interest in computational technologies that combine the flows of real data and hydrodynamic forecasts using mathematical models. This is especially true for 4D technologies-the combination of the flows of observational data and forecasts in a certain spatio-temporal domain. These methods have received the greatest applications in meteorology and oceanography, where observations are assimilated into numerical models. The purpose of assimilation procedures is to construct or refine the initial and boundary conditions (or other model parameters) to improve the accuracy of a prediction model [24][25][26][27][28][29][30][31][32][33][34].
At present, two main approaches are well known for the assimilation of observational data in models of geophysical hydrodynamics and oceanography. The first approach is based on the methods of probability theory and mathematical statistics. Historically, it is based on the classical least squares method proposed by C.F. Gauss (1795) and A. Legendre (1805); its rigorous justification and limits of applicability were given by A.A. Markov [35] and A.N. Kolmogorov [36]. From a methodological point of view, the least squares method can be considered the predecessor of the methods of optimal interpolation, the Kalman filter and their subsequent modifications, widely used in various fields of science and technology. The method belongs to the family of error theory methods used to estimate unknown quantities from measurement data, taking into account the random nature of measurement errors.
The second approach is based on the methods of calculus of variations, optimal control [28,37] and the theory of adjoint equations [4,27,[29][30][31]34]. Compared to the statistical method, the variational method has greater versatility. It allows, on a unified methodological basis, to solve the problems of initializing hydrophysical fields, assessing the sensitivity of a model solution, identifying model parameters, etc. The variational approach can be applied by assimilating information of various types and measuring systems. To solve these problems, the technique of 4D variational assimilation (4DVAR) is usually used. The main idea of the method is to minimize some functional that describes the deviation of the model solution from the observational data, and the minimum of this functional is sought on the model trajectories, in other words, in the subspace of model solutions.
The problem is formulated in a four-dimensional space-time domain and requires the solution of a coupled system of direct and adjoint equations in forward and backward time, respectively. Compared to the forecast Cauchy-type problem, the 4D data assimilation problem is much more complicated from the computational point of view. First, the problem adjoint to the original nonlinear problem has a more complex form. Secondly, when solving the adjoint problem, it is required to store a 4D array of the solution of the direct problem in machine memory. Third, when solving the problem of 4D data assimilation, iterative processes are used, and their convergence in many cases is under = question.
Ocean general circulation (OGC) models are very complex systems. Their basis is nonlinear differential equations describing the evolution of three-dimensional fields of currents, temperature, salinity, pressure and density [1][2][3][4]. In the operator of the ocean dynamics equations, two main parts can be distinguished. First, is an established classical basis-a subsystem that describes the dynamics of a rotating fluid in the framework of the approximations traditionally used in oceanology. Second comes the part which includes various kinds of physical parameterizations and is changed as our understanding of natural phenomena improves. In this regard, the formulation of the model and the development of an efficient numerical method are based on the splitting into physical processes. Using physical considerations, we select the main subsystems and then carry out their numerical solution [4,7,10].
The main feature of the INM RAS model INMOM, which distinguishes it from other well-known ocean models (see, for example, [1][2][3]), is that its numerical implementation is based on the method of splitting according to physical processes and spatial coordinates [10,38]. The complex operator of the differential problem of ocean dynamics is represented as a sum of simpler non-negative operators describing individual physical processes: three-dimensional transport-diffusion of moment, heat and salt, adaptation of density and current fields, dynamics of external gravitational waves, etc. This representation allows us to split the complex problem operator into a number of simpler ones and solve it in time using explicit or implicit schemes.
The aim of this paper is to present the technology for modeling and prediction of marine hydrophysical fields based on the 4D variational data assimilation technique developed at the Marchuk Institute of Numerical Mathematics, Russian Academy of Sciences (INM RAS). The technology is based on solving equations of marine hydrodynamics using multicomponent splitting, thereby solving an optimality system that includes adjoint equations and covariance matrices of observation errors. The hydrodynamic model is described by primitive equations in the sigma-coordinate system, which is solved by finite-difference methods. The technology includes original algorithms for solving the problems of variational data assimilation using modern iterative processes with a special choice of iterative parameters. The methods and technology are illustrated by the example of solving the problem of circulation of the Baltic Sea with the 4D variational data assimilation of sea surface temperature information.

Mathematical Model of Marine Circulation
A mathematical model of the sea dynamics is considered in geographical coordinates. Let U = (u, v, w) be the velocity vector, ζ the sea surface level function, T the temperature and S the salinity.
In the domain D of variables (x, y, z) for t ∈ (0,t) we consider the system of equations of hydrothermodynamics for the functions u, v, ζ, T, S under hydrostatics and Boussinesq approximations [39], with the Lame coefficients for a spherical coordinate system [40]: where is the forcing, f T , f S , f P are the functions of the "internal" sources, ρ 0 = const ≈ 1 is the mean density, T (0) , S (0) are the reference values of temperature and salinity, β TS (T, S) is the sum of all other terms of the expansion of the function of state ρ = ρ (T, S), f 3 ≡ f 3 (x, y, t) is the function related to the tide-generating forces, β T , β S , γ, g = const, A ϕ ϕ ≡ −Div â ϕ Gradϕ , Θ(z) ≡ (R − z)/R≈ 1 and R is the Earth radius. The first equation in (1) is a two-component vector equation for horizontal velocities u and v. The second equation describes the evolution of the sea surface level. The two last equations govern three-dimensional transport-diffusion process of heat and salt, respectively.
For system (1) in D × (0,t) the following boundary conditions [40] on the sea surface Γ S = Ω are set: where τ (a) x , τ (a) y , γ T , γ S , T a , S a , Q T , Q S are given functions, U (−) n = (|U n | − U n )/2, U n | z=0 = −w| z=0 . The first two equations in (2) define the momentum flux on the sea surface. The last two equations describe heat and salt flows through the sea surface. The vertical velocity w = w(u, v) is defined by the formula where r = R − z ≈ R. Equation (3) is obtained by integrating the continuity equation with respect to depth. In addition, boundary conditions are set on the "solid side wall" Γ w,c , on the "liquid part of the side wall" Γ w,op and at the bottom Γ H [33]. Thus, the normal component of the velocity vector is supposed to be zero on the solid side wall and at the bottom (the so-called "no normal flow" condition). The same is supposed for the normal components of the gradients ∇T and ∇S of temperature and salinity, respectively (no heat and salt fluxes on these boundaries).
Initial conditions for u, v, T, S, ζ have the form where u 0 , v 0 , T 0 , S 0 , ζ 0 are given functions.
The problem of large-scale sea dynamics in terms of the functions u, v, ζ, T, S is formulated as follows: find u, v, ζ, T, S, that satisfy (1)-(4). If the functions u, v, ζ, T, S are found, then the function w is determined by Equation (3).
Equations (1)-(4) are called primitive equations describing OGC. Primitive OGC models have several features. These are the Boussinesq and hydrostatics approximations (since the vertical size of the ocean is much smaller than the horizontal size, the equation for vertical velocity is replaced by the hydrostatics relation), the incompressibility of sea water with a free sea surface, the Coriolis force (the effect of the Earth's rotation) and the turbulent nature of the diffusion and dissipation process, taken into account. These features are reflected in the structure of governing equations; they require special physical parameterizations and searching for efficient numerical algorithms.
Note that the above boundary conditions may be modified depending on a specific physical problem. Problem (1)-(4) is approximated by the splitting method using the finite differences [10,38].

Splitting Method and Main Features of the Numerical Model
The main features of the numerical model of sea dynamics (1)-(4) are the simultaneous use of the splitting method [10,38] and the transition to the σ-coordinate system [7,10]. We used these two components in tandem to build efficient computer technology for 4DVAR ocean data assimilation.
Let us explain the connection between these two components. The essence of the splitting method is that the process of solving a complex problem is divided into several steps in which simpler problems are solved [38]. The splitting method is used to solve classical evolutionary problems, each equation involving a time derivative. However, the problems of geophysical hydrodynamics and ocean dynamics are described by non-classical systems of equations, some of which do not have time derivatives. Such equations are sometimes called equations of the diagnostic type. These include the hydrostatic equation-the third momentum equation and the continuity equation. The absence of several terms in the hydrostatic equation, including the time derivative, is due to the fact that the ocean domain is quasi-two-dimensional because the ocean depth is much less than its horizontal size [1]. The absence of a time derivative in the continuity equation is related to the incompressibility of sea water. Another difficulty of the ocean dynamics equations is the boundary condition on the sea surface. In modern models of the ocean general circulation it is formulated on a free moving surface that describes sea surface. Because of this, the ocean model domain varies in time. These difficulties can be circumvented by introducing the σ-coordinate system [7,10]. In this case, the ocean dynamics problem is formulated in the form of a system of evolutionary equations, and the problem solution domain does not depend on time: its horizontal boundaries do not change, and the vertical coordinate changes from zero to unity.
Let us introduce the grid on [0;t]: 0 = t 0 < t 1 < ... < t J−1 < t J =t, ∆t j = t j − t j−1 and consider problem (1)-(4) on (t j−1 , t j ), assuming that the vector of the approximate solution φ k ≡ (u k , v k , ξ k , T k , S k ), k = 1, 2, ..., j − 1 at the previous intervals, is already defined. To approximate the problem, we apply one of the schemes of the total approximation method [38], which consists of the implementation of the following steps (to simplify the notation, the index j for all components of the solutions of the subproblems is omitted at these steps).
Step 1. Consider the problem under corresponding boundary and initial conditions.
Step 2. Solve the problem under appropriate boundary and initial conditions.
Step 3. The system is solved under corresponding boundary conditions, and the function ξ j ≡ ξ (1) is taken as an approximation to ξ on (t j−1 , t j ). Then the following problem is solved: where is taken as an approximation to the exact vector u on D × (t j−1 , t j ), and the approximation w j ≡ w u j , v j to the vertical component of the velocity vector is calculated.
Thus, when steps 1-3 are implemented, after the first step we get an approximation to T, after the second an approximation to S, and after the third step we get an approximation to u = (u, v), ξ. Therefore, the subproblems at these steps are independent of each other and may be solved in parallel.
Note once again that for numerical solution of the complete problem (1)-(4) we use the σ-coordinate system. The transition to the σ-system can be carried out at the stage of considering the complete problem before applying suitable splitting schemes and other numerical procedures [41]. However, the transition to σ-systems is also possible after applying splitting schemes, i.e., in our case, as applied to problems (5)- (9). A number of other approaches for numerical solution of problems (5)-(9) are described, for example, in [39,41].

4D Variational Assimilation of Sea Surface Temperature Data
Variational data assimilation is based on the idea of minimization of a cost function associated with observational data on the trajectories (solutions) of the model under consideration. Thus, the assimilation problem is formulated as an optimal control problem in order to determine unknown model parameters [28][29][30][31][32][33][34].
Suppose that in problem (1)-(4) an additional unknown ("control") is the total heat flux function on Γ S . Let us introduce the cost function in the form where Q (0) = Q (0) (x, y, t) is a given function, α = const > 0 is the regularization parameter, T obs is the function of observations on the sea surface Ω and R is the observation error covariance operator. The problem of variational data assimilation is formulated as follows: find a solution to problem (1)-(4) and a function Q, such that functional (10) takes the minimum value.
The optimality system, which determines the solution of the formulated problem of variational data assimilation, includes the direct problem (1)-(4), the adjoint problem and an additional optimality condition of the form: The last equation means that the gradient of the functional J α (Q) with respect to Q equals zero. Here T * is the solution of the adjoint problem, which in the case of applying the splitting method is determined at step 1 in the form: Problem (12) is adjoint with respect to (5). As shown in [42], the optimality system that determines the solution of the formulated problem of variational data assimilation reduces to the sequential solution of some subproblems on t ∈ t j−1 , t j , j = 1, 2, . . . , J.
We formulate some of the algorithms for solving the problem under consideration. The construction of approximate solutions of the complete numerical model with the simultaneous determination of Q by variational assimilation of T obs may be carried out by the following iterative algorithm. If Q (k) is the already constructed approximation to Q on t j−1 , t j , then after solving the direct problem with Q ≡ Q (k) the corresponding adjoint problem is solved, and the following approximation Q (k+1) is computed by: with the parameter γ k , which is chosen so that the iterative process under consideration converges [42]. After computing Q (k+1) the solution of the direct and adjoint problems is repeated already with the new approximation Q (k+1) , and then Q (k+2) is calculated, and so on. Iterations are repeated until a suitable convergence criterion is met. According to [42], for γ k one can take the parameters which may significantly accelerate the convergence of the iterative process. The above-considered variational assimilation problem belongs to the class of four-dimensional variational data assimilation problems. The splitting method is considered here as a method of approximation of the original model, while the problem of variational assimilation itself is solved on the set D × (t j−1 , t j ), j = 1, 2, . . . , J.

Implementation
A detailed description of the numerical INMOM model can be found in [41] (see also [43,44]). To illustrate the operation of variational data assimilation algorithms, we performed numerical experiments for the problem of reconstructing the heat flux function Q in the Baltic Sea by variational assimilation of the observed sea surface temperature (SST). The calculations were carried out on the basis of the sigma-model of the Baltic Sea, developed at the INM RAS [13], based on the splitting method and supplemented by the variational assimilation of data taking into account covariance matrices of observation errors. The approximation of the model on the "grid C" [13] was used. The transport-diffusion equations were solved by splitting with respect to geometrical coordinates, using implicit schemes in time. For numerical experiments, modifications of the boundary conditions were chosen according to [10].
The model area of the Baltic Sea is located from 9.375 • E to 30.375 • E and from 53.625 • N to 65.9375 • N. The spatial resolution of the model is 1/16 • × 1/32 • × 25 in longitude, latitude and vertically. The grid area in the horizontal plane contains 336 × 394 nodes; the sigma levels are unevenly distributed in depth. The first point of the "grid C" [41] has the coordinates 9:406 • E and 53:64 • N. The mesh sizes in x and y are constant and equal to 0.0625 and 0.03125 degrees, respectively. The time step is 5 min. The computational domain with the topography of the Baltic Sea is presented in Figure 1. Here the section along 19 • E is shown (black dashed line), which was used in the experiments to present the depth values of temperature and salinity. The model was started from climatic temperature and salinity fields and zero initial velocities and ran with atmospheric forcing obtained from reanalysis until the solution reached a nearly stable, dynamically consistent regime, and after that the result of calculation was taken as an initial condition for further running of the model. The spinup time was about 1 year for the Baltic Sea. To start the 4DVAR assimilation procedure for the heat flux reconstruction, we consider the model forecast at the end of the previous time interval as an initial condition for calculations on the current time interval.
Sea surface temperature observation data provided by Copernicus Marine Service were used in the study. The daily mean sea surface temperature data of the Baltic Sea of the Danish Meteorological Institute, based on measurements of radiometers (AVHRR, AATSR and AMSRE) and spectroradiometers (SEVIRI and MODIS) [45], were used as T obs . To recalculate the observational data on the computational grid of the numerical model of the thermodynamics of the Baltic Sea, data interpolation algorithms were used [46,47].
The diagonal elements of the covariance matrix R, recalculated based on the statistical properties of the observational data, were taken as weighting coefficients in the cost functional when solving the data assimilation problem. Statistical characteristics (the mean and the variance) were computed based on observational data for 35 years, from 1982 to 2017, separately for each day of the year. Era-Interim global reanalysis data were used as a forcing. For Q (0) we used the mean climatic flux obtained according to the NCEP (National Center for Environmental Prediction) reanalysis. Note that Q (0) has a physical meaning here; it is not only an initial guess for Q but also a parameter calculated from atmospheric data and taken in the model for temperature boundary condition on the sea surface when the model runs without the assimilation procedure. Using the presented model of hydrothermodynamics, supplemented by the procedure of assimilation of the surface temperature T obs , calculations were performed in the Baltic Sea, in which the assimilation algorithm worked only at some time instants t j , with t j+1 = t j + ∆t j+1 . This means that until t j the calculation was done according to the model without the assimilation algorithm, and starting from t j the assimilation algorithm was turned on, and the initial assimilation field was set from the previous calculation for t = t j . When implementing the assimilation procedure at one time step (t j−1 , t j ) a system of the form (5)- (9) was considered, where by (5)-(9) we mean finite-dimensional analogues of the corresponding problems. The assimilation block diagram, according to the algorithm (13), is presented in Figure 2. It includes the data assimilation module based on minimization of the cost function (10), using the observation data T obs .

Results of Numerical Experiments for the Baltic Sea Circulation Model
Let us present the results of numerical experiments for the problem of reconstructing the heat flux function Q in the Baltic Sea by variational assimilation of the observed sea surface temperature data.
In Figure 3 the results of a numerical experiment are presented for a period of 15 days (January 2007). Figure 3a shows the mean value of the surface temperature field on 15 January according to satellite observations. The mean value of the sea surface temperature field on the fifteenth day of calculation (daily mean), obtained by the model without the assimilation procedure, is shown in Figure 3b. Figure 3c presents the mean SST field for the same day calculated with the use of the data assimilation procedure according to algorithm (13). From Figure 3a,b one can see that the model somewhat underestimates the surface temperature in the southern part of the Baltic Sea. This suggests that in the nearest future, attention should be paid to the parameterization of vertical exchange in the upper layer of the Baltic Sea. However, using the variational data assimilation procedure adjusts the model behavior and brings the calculation closer to observational data in the problem area (Figure 3c).    Figure 4 shows the differences in the daily mean SST values on the last day of the calculation (15 January). Figure 4a presents the deviation of the daily mean SST field calculated by the model without the assimilation block (T model ), from the observational data T obs . Figure 4b shows the deviation from the observation data of the SST field calculated with the use of the assimilation procedure (T assim ). The difference between the calculations using the model without assimilation and with the assimilation procedure, in which the elements of the covariance matrix of observation errors were used, is presented in Figure 4c. One can see from Figure 4 that the inclusion of the observational data assimilation block in the model significantly improves the prognostic properties of the model. It is interesting to note that the assimilation procedure did not affect the behavior of surface temperature in the northern part of the sea (Figure 4c). This part of the Baltic Sea is covered by sea ice in winter. To improve the accuracy of modeling, it is necessary to include functions describing the sea ice melt/formation in the data assimilation procedure of hydrophysical fields in this region.
(a)   Figure 5 shows the daily mean temperature on the 15th day in the vertical section of the Baltic Sea along 19 • E (this section is presented in Figure 1 by black dashed line). The differences between the model runs without assimilation and with assimilation in Figure 5c show that when assimilating SST, a change in temperature occurs only in the surface layer, and the temperature at the other levels remains almost unchanged. Thus, we can conclude that the depth values of temperature simulated by the model can be corrected only if in-depth observation data are available. We also note that the assimilation procedure, while adjusting the water temperature in the upper layers, makes it warmer when calculating in the winter period for the studied water area (Figure 5b). The results of numerical experiments also showed that the SST assimilation procedure has a weak effect on other components of the complete solution of the problem, such as salinity, current velocities and sea surface level. Figure 6 shows the daily mean salinity in the Baltic Sea on the 15th day of model integration. Figure 6a,b corresponds to a section along 19 • E. Here S model is the model salinity without the SST assimilation procedure, and S assim is the model salinity obtained using the assimilation of observational data. Note that the salinity in the calculation (see Figure 6a) generally corresponds to the hydrological regime of the Baltic Sea. Figure 6b demonstrates that the salinity remains almost unchanged after assimilation of sea surface temperature for 15 days. From a physical point of view, our problem is to estimate the heat flux on the sea surface using the synthesis of observational data and model data. It is clear that surface heat flux is mainly controlled by temperature of the upper layer and surface winds, so the observations of the SST should be taken into account first of all. Improving the simulation of other components of the marine circulation depends on the one hand on the improvement of physical parameterizations, and on the other hand on the data of observations of salinity, currents and sea level. As an additional control parameter, surface wind should be included with the SST.  Figure 7 shows the circulation in the Baltic Sea after 15th days of integration. For Figure 7a the calculation was carried out without the data assimilation procedure, and in Figure 7b the results of the experiment are presented for the same time period, but using the procedure of assimilation of sea surface temperature. The difference between these values is shown in Figure 7c. Note that the assimilation of only SST does not have a significant effect on the circulation of the Baltic Sea waters in short-term calculations, and the velocities differ by less than 5 percent. Note that the model adequately describes subsurface currents in the Baltic Sea (see Figure 7a), and the assimilation only slightly corrects them. According to Figure 7c, insignificant deviations of the velocity field are observed in a small subdomain of the southern part of the Baltic Sea in the vicinity of the the point 20 • E, 56 • N. Namely, in this region a significant change in sea surface temperature during assimilation is observed, which possibly affects the velocity.  Tables 1-3. Table 1 presents the differences in temperature, salinity and velocity values when calculating with and without procedure of variational assimilation of SST on standard depth levels for the point 20 • E, 63.2 • N in the Gulf of Bothnia. Similar results are given in Table 2 for the point 19.5 • E, 55 • N in the Gulf of Finland and in Table 3 for the point 24.3 • E, 60 • N in the southern part of the Baltic Sea. From the tables it follows that during variational assimilation for 15 days, the temperature difference can reach up to 1.5 degrees on the sea surface; however, it drops with depth. Thus, the assimilation of surface temperature most significantly affects the upper layers of the sea; the effect on the lower layers is mostly insignificant. The effect on salinity and water circulation during the assimilation of only SST is not significant in the case of short-term calculations.

Discussion
Let us recall that the main feature of the INMOM model is that its numerical implementation is based on the splitting according to physical processes and spatial coordinates. The splitting of equations can be multicomponent. The principal is splitting into physical processes: transfer-diffusion, adaptation and so on. Then, a secondary splitting of the main process can be carried out, for example, into the transport-diffusion of a substance along separate coordinates. Our experience shows that this method gives good results in modeling the dynamics of the seas and oceans [4,10,41].
The multicomponent splitting method allows for stable calculations, and it is convenient when writing codes for complex hierarchically developed systems. A natural property of a split model is its modular principle: a separate problem is related to a separate module. Moreover, when solving 4DVAR problems, the adjoint analogue is constructed for a separate subsystem, which simplifies the process of constructing an optimality system. Each subsystem can have its own adjoint. The total model can be "assembled" from a different number of modules. The computational characteristics of the model can be improved by changing individual computing modules. In the considered method of variational data assimilation, the model is used in the calculation at each iterative step of the optimality system. Therefore, at each calculated time instant, the solution of the variational assimilation problem satisfies the boundary and initial conditions of the model, and the numerical stability of the model is not violated. When running the model without assimilation with real forcing, it was noted that the model somewhat underestimates the temperature in the southern Baltic Sea in winter. Therefore, a variational assimilation method was proposed to correct the boundary condition for the possible improvement of the model properties. Using the algorithm proposed in the paper, it was possible to adjust the model behavior and bring the calculation closer to the observational data not only in the problem area, but throughout the water area. Note that the method with control at the boundary, when using the splitting method, requires less computational time-consuming compared to the more common method of variational initialization, since in this case the iterative process of solving the optimality system does not need to be carried out over the entire assimilation interval. We also note that it is possible to jointly use the initialization method and the method considered in the paper to improve the prognostic properties of the model.
Numerical experiments showed that the inclusion of the data assimilation procedure increases the calculation time by about 10%. This is due to the fact that the splitting method is used, which allows us to use implicit schemes and simplify the problem at each time step and thereby reduce the number of internal iterations of the assimilation algorithm. This gives the advantage to the variational approach over the statistical method, which requires ensemble calculations to construct covariance matrices at each time step.
The variational data assimilation approach has great versatility, because it is based on the fundamental principle to minimize a cost function describing the deviation of the model solution from the observational data, and the minimum of this functional is sought on the model trajectories. In the paper we have demonstrated the application of this approach to the problem of reconstruction of the heat fluxes, the boundary conditions on the sea surface. However, this approach allows one, on a unified methodological basis, to solve the problem of initialization of hydrophysical fields, to assess the sensitivity of the model solution, to identify model parameters, etc. Moreover, considering the weak formulation of 4DVAR can help to take into account the model errors and uncertainties via assimilation.

Conclusions
In this paper the methods and technology for modeling and analysis of marine hydrophysical fields based on 4D variational assimilation of observational data are presented. The developed methods can be used for monitoring and predicting sea currents, for solving problems of variational data assimilation in order to reconstruct boundary and/or initial conditions and other model parameters. The technology is based on a numerical model of marine circulation in a sigma-coordinate system, which includes a procedure of 4D variational data assimilation. Variational assimilation algorithms are based on modern iterative processes with a special choice of iterative parameters.
The main feature of the presented algorithms is the use of the splitting method for physical processes and geometric coordinates, which made it possible to simplify the problems under consideration at each time interval and make the implementation of 4D variational assimilation algorithms more efficient. Numerical experiments have shown the efficiency of the developed algorithms for 4D variational assimilation of sea surface temperature data using the example for reconstruction of heat fluxes on the sea surface. When developing and testing the algorithms and codes in this article, the Baltic Sea was chosen as a marine area for certainty. However, theoretical results and computational algorithms under consideration may be extended to other water areas.