The methods of data assimilation for sea dynamics problems are very important and are studied by many authors. We present our research on the development of variational data assimilation technique for sea dynamics problems with the aim to correct or to find the sea surface heat fluxes using the observational data. We consider the mathematical model of sea dynamics developed at the Marchuk Institute of Numerical Mathematics of the Russian Academy of Sciences (INM RAS). For observational data, we use daily mean SST observations from the Copernicus Marine Data Store. The variational data assimilation technique is based on a minimization of the cost function related to the observational data on the sea surface. The cost function involves the background and observation error covariance matrices. Our novelty is that the problem is reduced to a coupled system of model and adjoint equations, and it is solved by an iterative algorithm with the optimal iterative parameters. Numerical experiments for the Baltic Sea circulation model demonstrate the efficiency of the developed methodology. It is shown that the use of the developed data assimilation technique makes it possible to bring model calculations closer to observational data, and thereby helps to improve the forecast properties of the model.