1220 Ocean Dynamics (2019) 69:1217–1237 boundaries opened in the regions of main rivers. Further, HBM includes a sea ice model component that describes sea ice thermodynamics and incorporates Hibler-type dynamics (Hibler 1979). The BGC model ERGOM was originally developed by Neumann (2000) for the Baltic Sea and upgraded later by Maar et al. (2011) for the ecosystems in the North and Baltic Seas. ERGOM simulates the BGC cycling in the coastal seas using three phytoplankton groups (Cyanobacteria, Flagellates, Diatoms), two zooplankton size groups, four nutrient groups (nitrate, ammonium, phosphate, and silicate), two detritus groups (N-detritus and Si-detritus), oxygen and labile-dissolved organic nitrogen in the water column (lDON, Neumann et al. 2015). The phytoplankton and zooplankton groups are expressed in nitrogen concentrations. The chlorophyll a concentration and the Secchi depth are computed diagnostically (Doron et al. 2013; Neumann et al. 2015). Riverine load inflow of nutrients was derived from climatological data for major rivers. The boundary conditions for the BGC state variables are from the World Ocean Atlas (WOA05) as described by Maar et al. (2011). ERGOM is coupled one-way to HBM so that the physical fields influence the biogeochemistry, which itself does not influence the physics. 3 Data assimilation The data assimilation is performed using the ensemble- based error subspace transform Kalman filter (ESTKF Nerger et al. 2012b) provided by the parallel data assimilation framework (PDAF, Nerger et al. 2005, 2013), which are described in this section. 3.1 Parallel data assimilation framework The PDAF (Nerger et al. 2005, 2013 http://pdaf.awi.de) is an open-source software environment for ensemble data assimilation. It simplifies the implementation of the data assimilation system with existing numerical models by providing support to modify the model to compute ensemble forecasts and by providing fully implemented ensemble data assimilation methods. For the data assimilation, the model code is augmented by subroutine calls to PDAF. This changes the parallelisation of the model, so that it can simulate an ensemble of model states, which are then used in the analysis step of the data assimilation, where the observational information are incorporated into the model. 3.2 Error-Subspace Transform Kalman Filter The data assimilation method used here is the ESTKF (Nerger et al. 2012a). The ESTKF is an efficient variant of the ensemble Kalman filter, which uses an ensemble of Ne model states to represent the state estimate, as the ensemble mean, and its uncertainty by the ensemble spread. For an overview of different filter methods, see Vetra-Carvalho et al. (2018). The ESTKF performs a sequential assimilation by alternating forecast phases and analysis steps. In the forecast phase, all model states in the ensemble are integrated by the model until the time when observations become available. Then, the analysis step is computed in which the observational information is assimilated into the model states. Compared to the classical ensemble Kalman filter (EnKF Evensen 1994; Burgers et al. 1998), the analysis step of the ESTKF is a particularly efficient formulation because it takes into account that the number of the degrees of freedom for the analysis update is given by Ne ? 1, while the EnKF computes the update according to the usually much higher number of observations (see Nerger et al. 2005 for a comparison of the EnKF with the SEIK filter, which has the same efficiency as the ESTKF). Mathematically, the ensemble describes the degrees of freedom by spanning an error subspace of dimension Ne ? 1, which motivates the name of the filter method. In the analysis step, the ESTKF uses ensemble-sampled error covariances of the model forecast, the observation error and the observational values to estimate the true state of the system. The ESTKF does this as follows by computing transformation weights. Let Xk denote an ensemble matrix at time k in which each of the Ne columns represents one model state. The transformation of the forecast ensemble, Xfk into the analysis ensemble, X a k is given by the following: Xak = X?fk + Xfk Wk (1) where the overbar denotes the ensemble mean and Wk is a transformation matrix of size Ne × Ne. Given that the degrees of freedom given by the ensemble are Ne ? 1, this transformation matrix is calculated in an error subspace of dimension Ne ? 1 at time k. Below, we omit the time index k, as all calculations of the analysis step are at this time. The transformation matrix is computed as follows. First, the ensemble states are projected onto the error subspace by the following: L = Xf T, (2) where T is a projection matrix of size Ne × (Ne ? 1) given by the set of equations as follows: Tj,i = ? ???? ???? 1 ? 1 Ne 1 1? Ne +1 , for i = j , j < Ne ? 1 Ne 1 1? Ne +1 , for i = j , j < Ne ? 1? Ne , for j = Ne. (3)