A pre-operational three Dimensional variational data assimilation system in the North / Baltic Sea

This paper describes the implementation and evaluation of a pre-operational three dimensional variational (3DVAR) data assimilation system for the North/Baltic Sea. Univariate analysis for both temperature and salinity is applied in a 3DVAR scheme in which the horizontal component of the background error covariance is modeled by an isotropic recursive filter (IRF) and the vertical component is represented by dominant Empirical Orthogonal Functions (EOFs). Observations of temperature and salinity ( T/S) profiles in the North/Baltic Sea are assimilated in the year of 2005. Effect of the 3DVAR scheme is assessed by a comparison between data assimilation run and control run. The statistical analysis indicates that the model simulation is significantly improved with the 3DVAR scheme. On average, the root mean square errors (RMSE) of temperature and salinity are reduced by 0.2 C and 0.25 psu in the North/Baltic Sea. In addition, the bias of temperature and salinity is also decreased by 0.1 C and 0.2 psu, respectively. Starting from an analyzed initial state, one month simulation without assimilation is carried out with the aim of examining the persistence of the initial impact. It is shown that the assimilated initial state can impact the model simulation for nearly two weeks. The influence on salinity is more pronounced than temperature.


Introduction
In coastal and shelf seas, operational forecasting systems become feasible in recent years due to several reasons: increasing maturity of numerical models, advances in systematic and real-time monitoring, and progresses in data assimilation techniques and applications.At present, there are increasing Correspondence to: W. W. Fu (wfu@dmi.dk)requirements for operational forecasting systems in coastal and shelf seas.The objective of operaional forecast is to provide information in wide utility and availability to deal with marine environment problems and satisfy needs at different levels from the community (offshore oil industry operations, fish stock management, chemical contamination, and water quality associated with public heath, etc.).
In coastal and shelf seas, numerical forecasting is in itself a major challenge for the scientific community because of the specific and rich dynamics in the coastal regions.Though there have been some studies with various methods (Cossarini et al., 2009;Zhang et al., 2010;Barth et al., 2007;Kuropav et al., 2007;Dobricic and Pinardi, 2008;Fu et al., 2011), the coastal-shelf scale data assimilation is rather immature as compared to the global scale data assimilation and less applied in operational forecasting systems for a long time.This problem can be attributed to several factors.Firstly, one major problem is the lack of real time insitu observations in an operational sense; Secondly, quality of satellite observations such as SST and SSH is relatively poor in coastal waters than open seas because the data is more influenced by the cloud cover in the coastal regions.This renders the assimilation more dependent on in-situ observations; Thirdly, complex topography and coastlines also impose some technical constraints on coastal-shelf data assimilation schemes.In recent years, this situation has greatly changed with the rapidly increasing observations and knowledge about coastal processes.More observation are availbale with new coastal-shelf sea operational observation systems deployed.For example, observing systems have been largely improved in the European area.Among them, the Baltic Operational Oceanographic System (BOOS) is now providing real-time ocean observations and forecasts for the marine industry, public and other end-users.Its follow-up system InfoBoos provides online data delivery of both satellite and in-situ measurements over the Baltic sea since 1999.To some extent, the rapid developments of these operational S. Y. Zhuang et al.: A pre-operational 3-D variational data assimilation system observation systems make it possible to assimilate real-time ocean observations into a coastal-shelf forecasting system.In the second place, recent development of data assimilation techniques offers more choices to develop a coastal/shelf scale ocean data assimilation system.For the data assimilation applications in coastal-shelf forecasting systems, the highly inhomogeneous, non-stationary and anisotropic forecasting error covariance is one of the most important issues to be coped with.For that purpose, a complex data assimilation scheme such as 4DVAR may be preferable because the adjoint model can take into account the trajectory of background error evolution.However, such scheme requires prohibitive computational resources, and also demands huge manpower to develop and maintain the system due to its difficulty and complexity in technical implementation.As compared to the 4DVAR, 3DVAR or Optimal Interpolation (OI) is more computationally efficient.Regarding these schemes, how to deal with the high inhomogeneity in the coastaloffshore regions remains a serious issue.Nevertheless,due to the low cost of implementation and computation, as well as being able to use various types of observations globally, 3DVAR is still an appealing scheme for the coastal operational forecasting systems.In addition, when developing variational data assimilation systems at operational centers, 3DVAR has been a necessary prerequisite to the ultimate goal of 4DVAR assimilation algorithms.
In the North Sea and Baltic Sea, Danish Meteorological Institute (DMI) has been running a regional oceanographic physical model DMI-BSHcmod (Dick et al., 2001) and provide operational forecasts since 2001.After several years' operational practices and research activities, this coast/shelf forecasting system is required to be improved in several branches.Among them, one of the most pressing motivations is to develop a suitable data assimilation technique so that the real ocean observations can be utilized to improve the forecasting skill of the system.As for the North/Baltic Sea, the spatial and temporal resolution required to make realistic predictions are generally much higher than the resolution required for the adjacent deep ocean.Processes such as tides, breaking of internal waves and the barotropic response to high-frequency atmospheric forcing have time scales of hours and horizontal scales that can be of order 100 m or less.Moreover, the complex bathymetry and topography requires some special treatments in the assimilaion schemes in order to realistically reproduce the background error covariance, which is apprently not satisfying the assumption of being Guassian (Fu et al., 2011).A 3DVAR scheme using anisotropic recursive filter (ARF, Hayden and Purser, 1995;Purser et al., 2003a) is tested in (Liu et al., 2009).Its capability of modelling the anisotropic background error covariances structure helps to certain extent make use of the subsurface in-situ observations.Fu et al (2011) also apply an Ensemble Optimal Interpolation (EnOI) to assimilate temperature and salinity profiles.A simplified Kalman filter was implemented in the operational DMI-BSHcmod since 2006 (Larsen et al., 2006), but this scheme is heavily reduced to only address the SST data.To use temperature and salinity profiles, we need to either expand this method or implement a suitable scheme according to our practical need.In consideration of computational cost and technical maintenance, a 3DVAR scheme is chosen for a pre-operational data assimilation system.Both IRF and ARF are applied in the 3DVAR data assimilation system.But as a preliminary step, only IRF is employed to carry out the experiments for the system verification and validation in this study.ARF will be adopted and tested in the next step for future practical applications.In this implementation, design for scalability of observation operators is taken into account so that the increasingly expanded ocean measurements from various platforms can be more easily added into the data assimilation system.
The rest of the paper is arrangesd as follows: Section 2 describes the basic 3DVAR scheme; Section 3 presents the preconditioning and transform of variables in 3DVAR.Section 4 describes the operational forecasting model and observations used for assimilation experiments.Section 5 presents the experiment with a single observation.Results from 1 year assimilation experiment with T/S profiles are given in Section 6. Section 7 gives a conclusion and an outlook on the future work.

General formulation
In general, the basic scheme of 3DVAR is to find the optimal solution of the model state x which minimizes the following cost function (e.g., Lorenc, 1997;Dobricic and Pinardi, 2008): x is the model state to be estimated.It usually refers to analysis state vector.x b is the background state vector, y o is the observation state vector.H is the non-linear observational operator with which the analysis equivalent of observation y = H (x) can be obtained and compared with the observation measurements.The superscript T denotes matrix transpose.
In the cost function, the misfit between analysis and background is weighted by the background error covariance B, and the misfit between analysis and observation is weighted by the observational error covariance R. Usually the optimal solution is found by minimizing the cost function J (x) with respect to x, in which its gradient is also needed for determining the search direction and iteration steps in the minimizing algorithm: Ocean Sci., 7, 771-781, 2011 www.ocean-sci.net/7/771/2011/Following an incremental method (Courtier, 1994), Eq. ( 1) is linearized around the background state into the following form: where d = y o − H (x b )is the innovation vector, H is the linearized observation operator evaluated at x = x b and δx = x − x b is the analysis incremental vector.In this way, the original problem translates into finding an incremental analysis δx.Equation (2) becomes: In the current scheme, the state vector contains only temperature and salinity : x = [T S] T (5)

Numerical algorithm of minimization
For a typical coastal ocean data assimilation system, the order of the original size of background error covariance matrix B is about 10 6 ∼ 10 7 .After some applications of simplifying procedures like preconditioning variable transform and applying EOFs, the size of B is still very large, so minimization algorithm with high efficiency is crucial to solve the 3DVAR problem.
An available quasi-Newton L-BFGS algorithm (Byrd et al., 1995) is applied in this 3DVAR system to minimize the cost function (Eq.6).Indeed, L-BFGS uses a limited memory variation of the Broyden-Fletcher-Goldfarb-Shanno (BFGS) update to approximate the inverse Hessian matrix (denoted by H K ).Unlike the original BFGS method which stores a dense n × n approximation (if the lenghth of state vector is n), L-BFGS stores only a few vectors that represent the approximation implicitly.Due to its moderate memory requirement, L-BFGS method is particularly well suited for optimization problems with a large number of variables

Preconditioning and transform of control variables
In practice, the most difficult issue in minimizing the cost function of 3DVAR is how to deal with the inverse of the background error covariance matrix B due to its tremendous size.One practical solution to this problem is to perform a preconditioned control variable transform (defined by δx = Uv) in the process of minimization (e.g.Lorenc, 1997;Weaver et al., 2003).Normally, the transform U is chosen to approximately satisfy the relationship B = UU T and when the control variable vector v is chosen, errors of the variables are assumed to be uncorrelated.With this assumption, the cost function is rewritten as follows: Equation ( 4) then becomes: In this way, the minimization can be carried out without handling the inverse of B. Although the optimal solution can be found with Eqs. ( 6) and ( 7) theoretically, further simplification is necessary since U still has a large size.Suppose the background error covariance between any two points can be separated into a product of vertical component by horizontal component, it is easy to demonstrate (Lorenc, 2000) that the computation of B implicitly involves the transform of U which includes a sequence of linear operators: where U H is the horizontal part of the control variable transform related to the horizontal mode of B, U V is the vertical part of the control variable transform related to the vertical mode of B, and U P is the physical transform related to the multivariate dynamic or physical constraints (e.g. the relationship between sea surface height, SSH, error and temperature/salinity error; the relationship between current and SSH etc.).The formulation and meaning of each linear operator is described in the following subsections.

Horizontal part of control variables transform
The background error covariance in our scheme is assumed initially to be isotropic and homogeneous Gaussian type.
Usually the Gaussian type spatial correlation can be efficiently modeled in two ways.One way is to repeatedly apply a Laplacian operator which is also the solution of the horizontal diffusion equation.This method is mostly used in those 3DVAR data assimilation schemes applied in global NWP spectral models and also in global oceanographic data assimilation (e.g., Weaver and Courtier, 2001).The other way is to apply a recursive filter to approximate the Gaussian correlation function (Lorenc, 1992;Hayden and Purser, 1995;Dobricic and Pinardi, 2008).While the recursive filter is widely applied in the meteorological data assimilation, applications are less documented in the oceanic literature.
The recursive filter operator is defined as follows: where A i is the original value at grid i, B i is the value after filtering for i = 1 to I , C i is the initial value after one pass of filter in each direction from i = I to i = 1.α is the filter coefficient to be determined.According to Hayden and Purser (1988), the boundary condition for one filtering pass is: The filter coefficient can be deduced through comparing the spectral response function of the filter with the correlation function.For a typical Gaussian correlation function: where ε 2 is the variance, x is the grid spacing and L is the horizontal correlation length scale, we have: where E = N x 2 /L 2 in which N is the number of filter passes.It is important to note, the IRF should be multiplied by a factor β = ε 2 √ 2πL in the application.For those error correlations being closer, the modeling of the Gaussian type correlations with IRF is not adequate.An alternative method called the second-order auto-regressive model (SOAR) can be chosen as the correlation function: In this case, the correlation coefficient is also calculated using Eq. ( 12), but with different E and β : In our scheme, both correlation models are employed in the 3DVAR for IRF as options in which L = 35 km, N = 10.The standard deviation of background error ε is set to be 0.5 for both T/S fields.Because R f is a symmetric and selfadjoint operator, we can have: It means that U H can be calculated by total of N/2 filter passes.

Vertical part of the control variables transform
Like many implementations of 3DVAR (e.g.Bark et al ., 2004;Dobricic and Pinardi, 2008 ), to reduce the computation of the background error covariances, an empirical orthogonal functions (EOFs) decomposition of the vertical component of background error covariance is used for the vertical transform.The vertical transform is written in the following form: Where E is a matrix of eigenvectors, while is a diagonal matrix with eigenvalues of the EOFs.The vertical component of background error covariances can be estimated either by empirical functions or from statistics of large-size samples of model forecasts.But as the first step toward an operational assimilation system, the vertical correlation is only modelled with a revised empirical function in the form of: where L z = 20 m is the vertical correlation scale, z is the thickness between any two layers.The structure of the vertical correlation function is shown in Fig. 1.Since only the most dominant EOFs are used, the number of EOFs is much smaller than the size of the control vector.Thus, the computational effort is significantly reduced as the size of the background error covariance is greatly decreased.
4 The implementation

DMI-BSHcmod
The model used in this study is the DMI operational model DMI-BSHcmod, which is a hydrostatic three dimensional circulation model developed by German Bundesamt für Seeschifffahrt und Hydrographie (BSH).It is a free-surface with hydrostatic and Boussinesq approximation.The model solves the Navier-Stokes equations for the currents U/V, temperature T and salinity S. The finite difference method is adopted for its spatial discretization in which the staggered Arakawa C grid (Arakawa and Lamb, 1977) is applied on spherical coordinates horizontally, and z-coordinate is applied vertically.An upwind time integrating scheme is used for handling the advection of momentumwhile a conservative, fully explicit, multidimensional scheme is used for advection of the temperature and salinity.The model has been running operationally at BSH since mid-90s (Dick et al., 2001)  The model is forced by hourly meteorological forcing (10 m winds, 2 m air temperature, mean sea level pressure, surface humidity and cloud cover) based on the DMI operational weather model HIRLAM (High Resolution Limited Area Model).The forcing has a horizontal resolution of about 15 km.The surface heat flux is parameterized using bulk quantities of both the atmosphere and sea or sea ice, respectively.The evaporation flux is taken into account only in the heat budget.The changes of water volume due to evaporation, precipitation and ice formation are ignored.River runoff is the daily averaged data derived from river measurements for 5 German rivers, hydrological simulations for 42 Baltic rivers and climatology values for the rest of the rivers.
The turbulence vertical mixing scheme is based on a k −ω turbulence model extended for buoyancy affected geophysical flows (Umlauf et al., 2003) but with a new set of coefficients.Through parameterizations it takes into account breaking surface waves.The atmospheric forcing of the turbulence model is provided through a new set of surface flux boundary conditions for k and ω.Different algebraic stability functions are applied for the vertical diffusivities of momentum, heat, salt and passive tracers (Canuto et al., 2002) and these have been reconstructed into new, computationally sound expressions in turbulent time scale, temperature gradient, salinity gradient, resolved velocity shear and unresolved but parameterized shear from breaking internal waves.

Coastline treatment
When RF is applied in 3DVAR, its boundary condition (see Eq. 10) is usually defined in regular grid points.Since the existence of coastline in regional ocean model, how to deal with the boundary conditions of RF becomes a problem.Dorbricic and Pinardi (2008) handle this problem by inserting some imaginary water points at the land near coast and then the RF can be calculated on the extended grid.Liu et al. (2009) use a "mirror reflection" boundary condition for solving this problem.The purposes of both methods are to make the RF computation more efficient.However, because the initial value of the analysis increment is zero in the incremental 3DVAR, the RF algorithm can be extended to the land grid points when it is applied to the full grid.The compuational expense is reduced by this treament, but the accuracy may be degraded for the extreme cases such as a very thin peninsula.

Observations
T/S profile observations are gathered from several sources.The first database is the ICES (International Council for the Exploration of the Sea) which provides CTD profile data in the Baltic Sea.The second database is MADS (The National Database for Marine Data), which is supplied by NERI (National Environmental Research Institute, Denmark).the third source is from the BSH which provides T/S profiles from fixed buoy stations in the North Sea/Baltic Sea.The last data source is the ferry boxes in the national and regional monitoring programs.
Figure 2 shows the spatial and temporal distribution of T/S profiles in the North Sea/Baltic Sea in 2005.The upper panel shows the spatial coverage of the locations of stations or sites packed from different resources.The total number of observations sites is 4099.The spatial coverage of the observations indicates an irregular distribution with relatively denser data in the North Sea and Danish transition waters, but with sparse data in the Baltic Sea, especially in its northeast part.The bottom panel gives the daily number of T/S profiles throughout 2005 and displays quite uneven distribution of measurements with time.
Most of the T/S profiles have already gone through a primary data quality control.For further application in data assimilation, we have designed a simple intrinsic quality control scheme in the 3DVAR in order to prevent the potential "poor quality" data from destroying the analysis.As mentioned in Section 2.1, the innovation vector, denoted by the difference between the background field and the observations, is used as an indicator.If the innovation exceeds a certain number of observation error standard deviations, the observation is discarded.The criteria are set up empirically based on our past validation results of the model.For example, the observation is discarded if the magnitude of the innovation is large than 2.5 • C or 2.0 psu.
Figure 3 shows the observation usage in our experiments.Applying the data quality control scheme, out of total number of 82 354 temperature and 79 148 salinity measurements, about 92.63 % temperature and 91.16 % salinity measurements can be used into the data assimilation system.www.ocean-sci.net/7/771/2011/Ocean Sci., 7, 771-781, 2011

Observation error covariance
The observation errors at different sites are assumed to be uncorrelated in space, thus the observation error covariance matrix R is diagonal and the specification of R is then reduced to a list of observation error standard deviations for each data type, variable, and level, etc.In our current scheme, all T/S profile measurements are packed together, so only the unitary standard deviation for T/S are specified for our primary experiment, The standard deviation for T is set to 0.5 • C and the value for S is 0.5 psu.More realistic estimation of observation standard deviation values would be taken into account in the future.

Test with isolated observation
Single observation test is an important step to evaluate the performance of 3DVAR in a complex system.It is helpful in investigating the basic feature of the background error covariance and testing the validity of the DA system.The solution obtained from a complex system should be close to the analytic solution if the DA system is properly implemented and effective.Suppose one observation is located exactly at a model grid point, according to the characteristics of the analytic solution for 3DVAR, the increment at this point can be simply obtained (e.g., Zhuang et al., 2005) with the form: where δx is the univariate analysis increment at the grid point, dis the innovation corresponding to the single observation, σ 2 B and σ 2 o are background error variance and observation error variance, respectively.
A single observation analysis was carried out with the 3DVAR scheme and an isolated temperature observation (11.6 • E, 56 • N) at the depth of 5 m in Danish Waters.Figure 4 shows the horizontal and vertical distribution of the analysis increment field.In the horizontal, the numerical solution presents a Gaussian-shape structure.In the vertical, a Gaussian vertical structure is shown, which can be deduced from its corresponding correlation functions.The maximum value of the increment, from Fig. 4, is 0.52 • C. By Eq. ( 16), the analytic solution can be calculated at the observation point, which is about 0.5 • C. The two solutions are not identical because interpolations in the DA system have effect on the solution.Moreover, influence of iterations for minimization is not negligible.

Experiments
Three experiments are performed in order to evaluate the impact of the 3DVAR data assimilation scheme on the model simulations.The first two experiments are carried out for the year of 2005.The first experiment is a control run, in which the model perform daily 24-h forecasts without data assimilation.The second experiment is the same as the first one except that all collected observations are assimilated into the model by 3DVAR scheme.In order to examine how long the impact of the assimilated initial state can persist, the third experiment is run only for one month starting with an assimilated initial state from the second experiment on 1 March 2005.No observations are assimilated during this period.The assimilation is carried out daily provided that there are observations.

Experimental set-up
In the experiments, the model is set up with two-way nested domains in which the coarse grid covers the North Sea/Baltic Sea area with horizontal resolution of 10 in longitude and 6 in latitude (approximately 11 km by 11 km), and the fine grid covers the Danish Waters and German Bight with horizontal

Results
Normally, 3DVAR schemes can be assessed by comparing the statistics of differences between the model state and in situ observations for the experiments with or without data assimilation.In our system, the innovations are calculated and saved in log files, which provide an easy way to estimate the RMSE and bias from those differences.For the data assimilation experiment, it should be noted that the innovations are calculated before the assimilation analysis time and www.ocean-sci.net/7/771/2011/Ocean Sci., 7, 771-781, 2011 the corresponding observations are not used.Therefore, the evaluation of the analysis quality can be regarded as independent (Dobricic and Pinardi, 2008).Additionally, the RMSE and bias of analysis minus observation (analysis residual) are calculated in order to show to what extent the analyses agree with observations.Figure 7 shows the spatially averaged RMSE and bias estimated from the innovation and analysis residuals of both temperature and salinity for the period January 2005-December 2005.For both the temperature and salinity, the RMSE of innovation is generally smaller throughout the year of 2005 when the data assimilation scheme is applied, indicating that the data assimilation improves the simulations.Particularly, the RMSE and bias is significantly reduced with the 3DVAR data assimilation from the middle of February to the end of March.
One of the foundations for 3DVAR is the unbiased assumption.Ideally, the background, observations as well as the analysis are all unbiased.In other words, most data assimilation systems are bias-blind, which are designed to correct random errors only (Dee, 2005).In this case, we have where the point bracket represents spatial or temporal mean.However, in reality, a systematic bias inevitably exists in both the background and observation space, and consequently in the analysis space.Therefore, the mean of analysis increment stands for the analysis bias while the mean of innovation represents the model bias, if the observation is assumed to be the "truth".In Fig. 7b and d the daily model biases are shown for both experiments.The results indicate that the model biases of temperature and salinity are obviously reduced with 3DVAR data assimilation.Figure 8 presents the overall statistics of the RMSE and bias calculated from the 3 experiments.On average, with the 3DVAR data assimilation, the temperature RMSE is reduced by 0.22 • C, and the salinity RMSE is reduced by 0.41 psu.The T/S bias is reduced by 0.23 • C and 0.3 psu, respectively.From above estimations, it can be found that the temperature and salinity of the model have positive biases.Relatively, the bias of salinity is more decreased after the assimilation.
Further investigation has been done by analyzing the time averaged profiles for T/S RMSE and bias calculated from the above experiment results.Because of uneven depth of the ocean bathymetry in the North and Baltic Sea, as well as the scarcity of observations in deeper ocean, the calculation is performed only covering the water volume above 100 m.Most of the observations are located above 100 m, making the statistics more reliable.Figure 9 shows the comparison of RMSE and bias of T/S innovation for the control run, assimilation run and the analysis residual.It is obvious that the data assimilation significantly reduces the RMSE and bias of the model simulation throughout the layers.For temperature, the RMSE is largest at the depth of 60 m, which corresponds to the typical depth of thermocline.The RMSE of salinity shows a clear declining tendency except at the depth of 40 m, where a spike is found.On average, the RMSE of temperature and salinity is reduced by 0.2 • C and 0.25 psu.Aside from the RMSE, the 3DVAR is shown to be effective in reducing the model biases of temperature and salinity (right panel in Fig. 9).The positive bias of the temperature is found with two maxima locating at the depth of 15 m and 60 m, respectively.For salinity above 60 m, the bias is positive with values up to 0.8 psu.Below 60 m, however, the bias is negative.After assimilation, the bias of temperature and salinity is reduced by 0.1 • C and 0.2 psu on average.
In order to examine how long the influence of the assimilated initial state will persist, simulation of one month has  During this period, no observation is assimilated but with an initial state obtained from the above data assimilation experiment.In this experiment, only the initial state contains the previous observational information.In Fig. 10, the result of this experiment is compared with the control run and assimilation run.From Fig. 10a, it is seen that the RMSE of temperature is clearly smaller than the control run at the first 2 weeks.For salinity, the RMSE is reduced by 0.5 psu at the first 2 weeks with the initial state taken from the assimilation run.Different from the temperature, RMSE of salinity is still smaller than the control run from 20 March to the end of the month.This provides useful clues when the time window of observations is chosen for an ocean data assimilation system.A daily cycle of assimilation may not be necessary.A time window of about 1 week may be used by assembling observations.This is also important for operational data assimilation when real-time observations are not available for T/S profiles.
From Fig. 10, it can also be noted that the bias of temperature is nearly the same starting from different initial states.Comparatively, the bias of salinity witnesses clear reduction of about 0.5 psu.The salinity experiences less variations than the temperature in this short period, which helps to retain the influence of the initial state.The RMSE and bias is also affected by the starting date.More experiments are needed in the future to further examine this issue.

Conclusions and discussions
The development and implementation of a 3DVAR data assimilation system for operational regional model are described in this paper.Similar to some other 3DVAR implementations, the most common approaches such as the incremental method, IRF and vertical EOFs are applied in our implementation.A single observation test is firstly performed to validate the effectiveness of the 3DVAR scheme in a complex model system.In order to further assess the method, three experiments for the year of 2005 have been carried out.The results show that the T/S simulation is improved as the RMSE and bias of the model results are greatly reduced with the data assimilation.On average, the RMSE of temperature and salinity is reduced by 0.2 • C and 0.25 psu.In addition, the bias of temperature and salinity are also decreased by 0.1 • C and 0.2 psu, respectively.Moreover, it is demonstrated that the assimilated initial state has a positive impact on the model simulation.From our experiment, the impact can persit for nearly two weeks.Relatively, the impact on salinity is more pronounced than temperature for the averaged RMSE.This offers a possibility of properly using observations in delayed mode (opposite to real-time) in an operational data assimilation scheme in the future.It suggests that we can also expect a better model simulation by conducting data assimilation once every 2 weeks when there are no real-time observations.The preliminary results are encouraging in some aspects, but there are still some issues to be addressed in order to further improve this data assimilation system.For instance, the ARF is worthy of replacing the IRF for modeling the anisotropic structure of the background error covariances in the North/Baltic Sea though it is computationally more expensive.The estimation of the background error could be further investigated with long time model simulations.Moreover, the multivariate control variables analysis could be a better choice for producing more consistent analysis and forecast.Finally, surface observations such as the sea surface temperature (SST) and sea surface height (SSH) will play a complementary role in the assimilation when T/S profiles are assimilated.Assimilation of satellite data in the coastal region will be important and beneficial to improve the preoperational assimilation system.

Fig. 1 .
Fig. 1.Vertical correlations of the background error covariance.Each line labeled with different colors stands for a vertical correlation between one layer located at the right top of the line and others.
and at DMI since 2001.It has been significantly improved in recent years with the support from the European Framework Program (FP) projects.As the salty bottom water flowing into the Baltic Sea originates from the northeast (NE) Atlantic, the model covers a large area (the Baltic/North Sea and NE Atlantic) in order to study Baltic/North Sea water exchange through the Danish Straits.Three nesting levels are applied: a 2-dimensional (2-D) NE Atlantic model (NOAMOD) is used to provide surge boundary conditions for a Baltic/North Sea 3-D model (NBCMOD).A 3-D coastal model covering the inner Danish waters (KUCMOD) is twoway nested with the NBCMOD (Fig. 1).All the three models are set up horizontally in spherical coordinates and vertically in z coordinate.Both the NOAMOD and NBCMOD have a horizontal resolution of about 6 nautical miles (nm) while the fine grid model has a horizontal resolution of about 1 nm.The 3-D models have in total 50 vertical layers.The top layer thickness is selected at 8 m in order to avoid tidal drying of the first layer in the English Strait.The rest of the layers in the upper 80 m have 2 m vertical resolution.The layer thickness below 80 m increases gradually from 4 m to 50 m.

Fig. 2 .Fig. 3 .
Fig. 2. T/S profile distribution in the North Sea/Baltic Sea during 2005.The upper panel shows the spatial coverage.The outer rectangle frame covers the model coarse grid domain and inner red rectangle frame covers the model fine grid domain.The bottom panel shows the temporal distribution.The blue column indicates the model coarse grid domain, the red column for the model fine grid domain.

Fig. 4 .Fig. 5 .
Fig. 4. The analysis increment corresponding to 1.0 unit innovation with an isolated T at Danish Water.Upper panel: horizontal distribution at 5 m; bottom panel: cross section at 56 • N.

Fig. 6 .
Fig. 6.DMI-BSHcmod bathymetry and domain setup.The bathymetry shown in the figure is for the coarse grid.The red rectangle frame outlines the inner fine grid.

Fig. 7 .
Fig. 7.The daily spatially averaged RMSE and bias of innovations (background minus observation) for the experiment with 3DVAR data assimilation (blue line), and the experiment without data assimilation (red line).The corresponding rmse and bias of analysis residual (analysis minus observation) are shown by green line.

Fig. 8 .
Fig. 8.The gross RMSE (left panel) and bias (right panel) of T/S for control (T ct /S ct ), 3DVAR data assimilation innovation (T in /S in ), as well as 3DVAR data assimilation analysis residual (T an /S an ).

Fig. 9 .
Fig. 9.The comparison of rmse and bias profiles for both temperature (A, B) and salinity (C, D) among the innovation from control experiment (red line), the innovation from 3DVAR assimilation experiment (blue line), and the analysis residual from 3DVAR assimilation experiment(green line).The experiments run daily forecast through the year of 2005.

Fig. 10 .
Fig. 10.The daily spatially averaged rmse and bias of innovation of both temperature (A, B) and salinity (C, D) for the experiments with 3DVAR data assimilation (blue line), without data assimilation (red line), as well as one month non-data assimilation run with assimilated initial value (black line) throughout March 2005.