Controlling atmospheric forcing parameters of global ocean models : sequential assimilation of sea surface Mercator-Ocean reanalysis data

Controlling atmospheric forcing parameters of global ocean models: sequential assimilation of sea surface Mercator-Ocean reanalysis data C. Skandrani, J.-M. Brankart, N. Ferry, J. Verron, P. Brasseur, and B. Barnier LEGI/CNRS, Grenoble, France Mercator-Ocean, Toulouse, France Received: 26 May 2009 – Accepted: 6 June 2009 – Published: 18 June 2009 Correspondence to: J.-M. Brankart (jean-michel.brankart@hmg.inpg.fr) Published by Copernicus Publications on behalf of the European Geosciences Union.


Introduction
One of the most accurate and ubiquitous information about the surface state of the ocean is provided by the satellite measurements of sea surface temperature.It is in particular significantly more accurate than the sea surface temperature that is simulated by any state-of-the-art general circulation ocean model.Part of this discrepancy is explained by the relative inaccuracy of the atmospheric parameters that are used to compute the air-sea momentum, heat and fresh water fluxes which determine the surface boundary condition of the ocean model (WGASF, 2000).There is thus an important potential benefit to expect from the improvement of these parameters using the available sea surface observations.In practice, the atmospheric parameters controlling the air-sea fluxes (i.e.air temperature, relative humidity, cloud fraction, precipitation or wind speed) are derived from atmospheric reanalysis data (as delivered for instance by the ECMWF or NCEP centers) and from a variety of satellite products.For instance, the atmospherically forced ocean hindcast simulations performed by The DRAKKAR Group (2007) compute their airsea fluxes by using forcing data that merge a variety of different data sets (in situ, satellite and NWP products), with objective corrections based on observations (Large and Yeager., 2008;Brodeau et al., 2009).Hence, as long as forced C. Skandrani et al.: Controlling atmospheric forcing parameters of global ocean models models are used to simulate the ocean component alone, the control of the atmospheric parameters using ocean surface observations is certainly an appropriate way of improving the realism of model interannual simulations, the accuracy of ocean reanalyses or the usefulness of sea surface temperature operational forecasts.It is thus also an important contribution to the GODAE1 objectives (GODAE, 2008), which is the reason why a large part of the MERSEA2 effort in the development of data assimilation has been devoted to this problem.
In this study, which has been conducted as part of the MERSEA project, this problem is investigated using idealized experiments in which Mercator-Ocean3 ocean reanalysis data are used as the reference simulation (i.e. the "truth" of the problem).Synthetic observation datasets (for sea surface temperature and sea surface salinity) are extracted from the reanalysis to be assimilated in a coarse resolution global ocean model.With respect to Skachko et al. (2009), who investigated a similar problem using twin assimilation experiments, the present study is thus more realistic, since the difference between model and reanalysis is now very similar in nature to the real error.It is closer to the real problem even if the experiments are still somewhat ideal in the sense that no real observations are assimilated, and that the full reference model state (the reanalysis, in three dimensions) is available for validation.Another difference with respect to Skachko et al. (2009) is that, in this paper, we extend the control vector to 6 atmospheric parameters instead of 2 turbulent exchange coefficients in their example (but we exclusively focus on the control of the parameters, while they also considered the joint optimal estimate of the ocean state vector together with the atmospheric parameters).However, in order to solve this more realistic problem, we needed to further develop the methodology towards a better specification of the prior information about the parameters and their associated uncertainty.We observe indeed that making appropriate assumptions on that respect is increasingly important as the estimation problem is becoming more realistic, because it is more and more difficult to make the distinction between forcing errors and the other potential sources of error in the system.An additional important objective is thus to find means of identifying properly the part of the observational misfit that can be interpreted as resulting from inaccurate atmospheric parameters.
In order to reach this objective, the plan is to apply sequentially a Bayesian inference method to compute piecewise constant optimal parameter corrections.A possible algorithm to solve this problem is to compute the optimal parameters by direct maximization of the posterior probability distribution for the parameters, using for instance a 4DVAR scheme (as done in Roquet et al., 1993or Stammer et al., 2004).But, in addition to the technical difficulties that the algorithm may involve, this solution requires that the cost function resulting from the optimal probabilistic criterion be quadratic or at least differentiable everywhere in parameter space, so that it is by no way straightforward to optimally impose strict inequality constraints to the parameters (by setting zero prior probability in prohibited region of the parameter space for instance).This is why, in this study, we prefer using a Monte Carlo algorithm to simulate the ocean response to parameter uncertainty, and use the resulting ensemble representation of the prior probability distribution to infer optimal parameter corrections from the ocean surface observations.It is in the specification of this prior probability distribution that two methodological improvements are introduced with respect to Skachko et al. (2009).First, the error statistics are computed locally in time for each assimilation cycle, by performing a sequence of ensemble forecasts around the current state of the system (while they are assumed constant in their study).And second, the probability distribution is assumed to be a truncated Gaussian distribution (as proposed by Lauvernet et al., 2009, as an improvement to the classical Gaussian hypothesis), in order to avoid the most extreme and non-physical parameter corrections.These two improvements are indeed found necessary to solve the more realistic assimilation problem at stake in this paper.
However, before explaining this in more detail, we first summarize in Sect. 2 the background existing elements that are used to perform the study: the ocean model, the assimilation method for parameter estimation and the Mercator-Ocean reanalysis data.Then, in Sect.3, we present the details of the method that is used to perform the assimilation experiments: experimental setup and statistical parameterization.And finally, in Sect.4, we discuss and interpret the results, focusing on the accuracy of the mixed layer thermohaline characteristics and on the relevance of the parameter estimates.

Background
In this section, we present the three existing ingredients that are used later as a background information to set up our assimilation system (Sect.3) and to perform the experiments (Sect.4): (i) the ocean model, focusing on the role of the atmospheric forcing parameters, (ii) the assimilation method, in order to introduce the various approximations and parameterizations that are needed to solve the problem, and (iii) the Mercator-Ocean reanalysis, from which the synthetic observations are extracted.

Ocean model
The OGCM used in this study is a global ocean configuration (ORCA2) of the NEMO-OPA model (Madec et al., 1998), using a 2 • ×2 • ORCA type horizontal grid, with a meridional grid spacing reduced to 1/2 • in the tropical regions in order to improve the representation of the equatorial dynamics.This is a free surface configuration based on the resolution of primitive equations, with a z-coordinate vertical discretization.There are 31 levels along the vertical, and the vertical resolution varies from 10 m in the first 120 m to 500 m at the bottom.The lateral mixing for active tracers (temperature and salinity) is parameterized along isopycnal surfaces, and the model uses a turbulent kinetic energy (TKE) closure scheme to evaluate the vertical mixing of momentum and tracers (Blanke and Delecluse, 1993).
The model is forced at the surface boundary with heat, freshwater and momentum fluxes.The fluxes through the ocean surface are estimated from the atmospheric parameters at the anemometric height, using the bulk semi-empirical aerodynamic formulas.Daily atmospheric variables (wind, humidity, air temperature, cloud coverage) from NCEP and monthly mean precipitation from CMAP (CPC Merged Analysis of Precipitation) are used to interactively diagnose the net heat and fresh water fluxes (Q NET and F W NET ), which can be written respectively: where Q S is the sensible heat flux, Q L the latent heat flux, Q LW the long wave radiation flux, Q SW the short wave solar radiation flux, and E, P , R are the three terms related to the fresh water budget, respectively evaporation, precipitations and river runoffs.The flux parameters which are involved in the computation of these quantities are the latent heat flux coefficient (C E ), the sensible heat flux coefficient (C H ), sea surface temperature (T w ), air temperature (T a ), air pressure, atmospheric specific humidity (q a ), wind speed (W 10 ), cloud coverage (C) and precipitation (P ).For more detail on these bulk formulas, the reader can refer to the CLIO (Coupled Large-scale Ice Ocean) model description in Goosse et al. (1999).
The turbulent latent and sensible heat fluxes are calculated from the classical ocean-atmosphere transfer equations (Large and Pond, 1982): -the latent heat flux: where ρ a is the air density, L e the vaporization latent heat, and q s the saturation specific humidity; -the evaporation fresh water flux : -the sensible heat flux: where c a p is the air specific heat; -the long-wave radiation flux, which is parameterized by following Berliand and Berliand (1952): where e a (in mb) is the vapor pressure deduced from q a , , the surface emissivity, σ sb the Stephan-Boltzmann constant, (1−χ C 2 ), a correction factor to take into account the effect of clouds; -the short-wave radiation flux, following the proposed formula by Zillmann (1972): where α is the ocean albedo, β, the zenith angle at noon and Q CLEAR , the solar radiation at the ocean surface in clear weather.
For the momentum flux, we did not use aerodynamic bulk formulas to calculate the wind stress vector.It is directly specified in the model, using ERS scatterometer data complemented by in-situ observations of TAO derived stresses (Menkes et al., 1998).No relaxation to observed SST and SSS is applied in our simulations.

Assimilation method
The purpose of this section is to briefly describe the assimilation methods that are applied to perform this study.Only general algorithms and equations are given here; the specific parameterizations on which they depend are presented in Sect.3.

Estimation of model parameters
The problem of estimating model parameters from ocean observations can be formulated using the Bayesian inference framework.From a prior probability p(α) for a vector of uncertain parameters α, and the conditional probability distribution p(y|α) for obtaining a vector of observations y given the vector of parameters α, the Bayes theorem: provides the posterior probability p(α|y) for the parameters given the observations.A best estimate α * for the vector of parameters can then be obtained as the mean (minimum variance estimator) or the mode (maximum probability estimator) of this posterior distribution.The most common methods to compute α * are direct minimization techniques (to compute the mode), Monte Carlo integration (to compute the mean) or a direct formula (for instance, if the distributions are assumed Gaussian).
In this problem, the observations y are usually not directly related to the parameters α, but to the model solution x that is a function of α, so that the probability distribution p(y|α) is usually defined as a function of the misfit between the observations and the model solution corresponding to α: y−Hx(α) (innovation vector), where H is the observation operator.This makes the computation of α * more difficult, either with direct minimization techniques, because every evaluation of the function to minimize requires one model simulation (and also one adjoint model simulation if the gradient is also computed), or with Monte Carlo methods, because they require an ensemble model forecast using an ensemble of parameter vectors drawn from their prior probability distribution.
In this study, Monte Carlo simulations are performed to compute the model counterpart x to an ensemble of parameter vectors, sampled from p(α).This ensemble forecast characterizes the prior probability distribution p(x) for the augmented vector x=[α, x(α)], characterizing the model response to parameter uncertainty.It is important to note that, up to this point, linearity has not been assumed, and that it is only at this stage that a Gaussian parameterization is used for the prior distribution p(x), with the consequence of linearizing the inference rules relating the model parameters α to the observations y (see below).In order to mitigate the effect of this linearization, the problem is divided in a sequence of short periods of time (assimilation cycles) and the parameters α are estimated separately and sequentially for every element of the sequence (see Sect. 3 for more detail).

Optimal estimate under Gaussian assumption
If the probability distribution p(x) and p(y|x) can be assumed Gaussian: where xb is the background simulation, P is the background error covariance matrix in the augmented space, Ĥ=[0, H] is the augmented observation operator and R the observation error covariance matrix, it is known that the posterior probability distribution p(x|y) is also Gaussian: where the mean xa and the covariance Pa are given by the standard linear observational update formulas: Equation ( 11) are also the equations of the observational update of a Kalman filter written for an augmented control vector (including model parameters in addition to the model state).This method can be applied to control other sources of error in addition to parameter error (as explained, for instance, in Skachko et al., 2009).However, in the present study, Eq. ( 11) are only going to be used to obtain improved parameter estimates.
It is interesting to note that the solution given by Eq. ( 11) is not equivalent to minimizing (where P α is the block of P corresponding to the vector of parameters) using a variational method, as soon as the function x(α) relating the model solution to the parameters is nonlinear.This variational solution only assumes Gaussianity of p(α) and p(y|x) while keeping the nonlinear function x(α) in the expression of which is not Gaussian.However, there is no prerequisite of Gaussianity in Monte Carlo methods, that may also offer other advantages.It is for instance easier to apply strict inequality constraints (by modifying the prior Gaussian assumption).This possibility is exploited in this study to confine the parameter estimates in a predefined region of the parameter space (see Sect. 3.5).

Reduced rank approximation
If the background error covariance matrix is available in square root form P= Ŝ ŜT , with a rank given by the number r of independent columns in Ŝ (the error modes), then the problem can be simplified to the estimation of a reduced vector ξ (of size r), giving the amplitudes of the correction to xb along each column of Ŝ, using a reduced observation vector η (of size r) resulting from the projection of the innovations onto the error modes Ŝ: where U (unitary matrix) and (diagonal matrix) are the matrices with eigenvectors and inverse eigenvalues of the r×r matrix: By transformation ( 14), the probability distributions (Eq.9) transforms to so that With the transformation U, the observational updates for every components of the ξ vector are independent (all matrices in Eqs.16 and 17 are diagonal).This is a simple way of obtaining directly the equation of the observational update for the SEEK filter (the reduced order Kalman filter developed by Pham et al., 1998), that are otherwise deduced from Eq. ( 11) using the Sherman-Morrison-Woodbury formula.

Mercator-Ocean reanalysis
Mercator-Ocean is an operational oceanography center based in Toulouse, France.It develops and runs operational ocean analysis/forecast systems specially designed to provide useful products for several downstream applications: research, institutional and operational applications, private sector applications and environmental policy makers.Mercator-Ocean also periodically delivers ocean reanalyses, that are produced using up-to-date ocean models, observations and assimilation methods.In this study, we are using the data from a Mercator-Ocean coarse resolution global reanalysis (PSY2G2).The main application of this reanalysis is to provide ocean initial conditions for coupled seasonal prediction applications (Balmaseda et al., 2008), but it is also used for research purposes as it provides a long coherent time series of the ocean state (from 1980 to present).The ocean model used in the reanalysis is similar in many points to the one that is used in our experiments: same numerical code OPA, same grid, same physics, but different atmospheric forcing, which is computed using the ERA-40 reanalysis for the period January 1979 to December 2001.The assimilated observations are subsurface temperature and salinity, SLA data and SST maps.The subsurface data come from the EN-ACT/ENSEMBLES data base provided by the CORIOLIS4 data center; the altimetric data are along-track SLA (from November 1992 to present) provided by SSALTO/DUACS; and the SST maps are obtained from the Reynolds OIv2 product (Reynolds et al., 2002).The reanalysis data assimilation scheme is a reduced order Kalman filter using the SEEK formulation (Pham et al., 1998).The forecast error covariance is based on the statistics of a collection of 3-D ocean state anomalies (typically a few hundred) and is seasonally variable.The statistical analysis produces temperature and salinity as well as barotropic velocity increments, from which zonal and meridional velocity fields are deduced using physical balance operators.For the present study, we extracted 1993 and 1994 data (temperature and salinity) for use in our assimilation experiments.

Setup of the assimilation experiments
The general idea of the experiments presented in this paper is to use the Mercator-Ocean reanalysis as reference simulation from which synthetic observation datasets are extracted, and to assimilate these observations into our ocean model as a constraint to the atmospheric forcing function.These experiments are ideal in the sense that no real observations are assimilated and that the full reference model state (in three dimensions) is available for validation.But they are also realistic (clearly distinct from twin experiments) because dif-ferences of model simulations with respect to the reanalysis are similar in nature to differences with respect to the real world.It is indeed expected that the assimilation of the real observations to build the reanalysis moved the model trajectory towards a more realistic description of the ocean.The starting date of the experiments is 30 December 1992 (with an initial condition from a standard simulation performed by Castruccio et al., 2008).The first six months are used as an initialization period for the assimilation system, so that the one year diagnostic period extends from 30 June 1993 to 29 June 1994.Figure 1 (left panel, top black line) shows the time evolution of the RMS error (difference with respect to the reanalysis) in the free simulation (i.e.without parameter corrections) for sea surface temperature (SST) and sea surface salinity (SSS), as computed over the world ocean south of 70 • N, to avoid some problematic ice covered regions.We can observe on the figure that the SST RMS error is stable in time (in the interval 0.85-1.05• C) but that the SSS RMS error is drifting from the beginning of the simulation (at a rate of about 0.1 psu per year).
The error is however very inhomogeneous horizontally, as can be seen in Fig. 2, showing maps of SST and SSS systematic error (top panels) and maps of SST and SSS error standard deviations (bottom panels), averaged over the one year diagnostic period.The largest systematic errors (up to 2 • C and 0.8 psu) or error standard deviation (up to 1.5 • C and 1 psu) are localized in the regions of the Western boundary currents and the Antarctic Circumpolar Current (ACC).These large errors are due to the poor representation and localization of the ACC and the boundary currents in our low resolution ocean model.Since these currents are associated to the most intense SST and SSS fronts in the ocean, it is mostly the misplacement of the currents that leads to the largest SST and SSS errors.In these regions, the atmospheric forcing function is not the dominant cause of error, so that the identification of forcing errors is almost impossible there (with a low resolution model).This is why they will be masked in several diagnostics involving horizontal averages.The right panel of Fig. 1, for instance, shows the same result as the left panel as obtained by masking the 10% of the ocean with the largest free run RMS misfit.The mask is thus different for SST and SSS, and can be visualized in Fig. 2, as the region with the largest misfit (i.e.essentially the Western boundary currents and a part of the ACC).The RMS error for this 90% subdomain remains quite large in the free simulation with a reduction of about 10% for SST and 30% for SSS with respect to the global RMS error.
The purpose of our experiments is to identify the part of the misfit between model simulation and reanalysis that can be explained by errors in the atmospheric forcing function.This means that the experiments are dedicated to the improvement of the forcing function of long term free model simulations: no correction is applied directly to the evolution of the model state, which remains a solution of the model dynamical equations.But we need to apply corrections on the 1993 1994 forcing parameters using a sequential assimilation method.The simulation is thus divided in a sequence of time intervals (the assimilation cycles, with a length of 7 days), for which we estimate the forcing parameters by combining optimally forcing prior knowledge and available ocean observations.In our experiments, the forcing parameter estimates are obtained using SST and SSS observations extracted from the reanalysis at the model resolution with global coverage and perfect accuracy.

Initial condition
In such a kind of experiment, an important difficulty occurs if we start the assimilation experiment from the initial condition of the free model simulation on 1 January 1993, because at the end of the first one-week assimilation cycle, most of the error (difference with respect to reanalysis) is due to initial condition error and not to errors in the atmospheric forcing (corresponding to this first cycle).To avoid this difficulty, we first present results that are obtained without initial condition error.In that way, a better view of the behaviour of the method can be produced, thus facilitating the interpretation of the results.The influence of initial condition error is only briefly considered in a second stage.
In order to build a perfect initial condition for the assimilation experiments, we simply perform an ideal assimilation experiment with perfect increment δx k , computed as the difference between the model forecast of the current cycle (number k) and the corresponding ocean state in the reanalysis.This increment is then introduced into the model using the incremental analysis update algorithm (as in Ourmières et al., 2006).Moreover, in order to distribute the effort over the first p assimilation cycles, this increment is divided by the factor max(p−k+1, 1).In that way, only one pth of the full increment is applied during the first cycle, and the cycle p is the first cycle with the full perfect increment.Figure 1 (left panel, dashed black line) shows the SST and SSS RMS error reduction during this initialization procedure from 30 December 1992 to 10 March 1993 (10 cycles of 7 days with p=8).This last date is the initial condition of our simplified assimilation experiments (with negligible initial condition error).For comparison purpose, a free model simulation starting from this perfect initial condition is also presented in Fig. 1 (lower black solid line), showing that the corresponding SST and SSS misfits quickly increase with time to reach asymptotically the typical misfit (for SST) or the typical trend (for SSS) that is observed in our original free model simulation (with wrong initial conditions).

Forcing parameters prior probability distribution
As explained in Sect.2.2, the estimation of model parameters using Bayesian inference requires the definition of a prior probability distribution for the parameters.And the first thing to decide about this probability distribution is the list of uncertain parameters to include in the control vector (step a in Table 1); other parameters are then considered perfectly accurate.In this study, we decide to estimate the following parameters of the atmospheric forcing function: air temperature (T a ), air relative humidity (q a ), cloud fraction (C), precipitation (P ), the latent heat flux coefficient (C E ) and the sensible heat flux coefficient (C H ). The reason for this choice is that they are expected to be the most important inaccurate parameters that are involved in the computation of the net heat and fresh water fluxes at the air-sea interface.They are assumed to be responsible for most of the error in the computation of these net fluxes, and thus to be one of the most important source of error in the heat and salt budget of the ocean mixed layer.One important missing parameter is wind velocity, which is a key parameter to control the heat flux computation (Mourre et al., 2008).The reason for which it is not included in the list of control parameters is that it is also correlated to the momentum flux (zonal and meridional wind stress components), which we have chosen not to control in this study.The decision not to control this important source of model error results from the necessity of proceeding step by step to avoid unpredictable difficulties in the solution of the inverse problem.However, we must be aware that, as any uncontrolled error (like the localization of the boundary currents), this can introduce compensation problems in the parameter estimates (see Sect. 4.3).
In order to define the prior probability distribution for these control parameters (step b in Table 1), we first assume that the error on the parameters is constant over the current assimilation cycle, which already means that the overal flux correction in our experiments is necessarily piecewise constant (with weekly forcing parameter increments).Second, we assume that the parameter error pdf is Gaussian N (0, P α ), with zero mean and with the covariance P α of the time variability of the parameters in the free model simulation (here over the period 1992-1998).Using time variability as a uniform way of parameterizing parameter uncertainties is obvioulsy rather crude, especially if a better information is potentially available (as for C E and C H ), but this is a useful simplification for the definition of statistics and the interpretation of the results.With this method, it is also easy to obtain directly a reduced rank parameterization of the covariance matrix.In the definition of P α , we indeed retain the 200 first EOFs of the full signal, representing about 92% of the total variance.Figures 3 and 4 show respectively the resulting mean and standard deviation maps for every parameter.As SST and SSS misfit standard deviations shown in Fig. 2, the parameters standard deviations are very inhomogeneous horizontally, and the patterns of maximum variability are very diverse.These figures are used as a reference in further discussions.
However, it must be mentioned here that most parameters are constrained by bounds: www.ocean-sci.net/5/403/2009/Ocean Sci., 5, 403-419, 2009 Table 1.The four steps of the parameter estimation scheme, with a short description of the procedure and basic assumptions.

Steps Procedure Assumptions
Step a: definition of the augmented control vector Include a list of uncertain atmospheric forcing parameters in the control vector: The other parameters are perfectly accurate.
Step b: forcing parameter prior probability distribution Postulate a Gaussian pdf for the atmospheric forcing parameters, based on their natural variability simulated by the free model over 7 years, 92-98.(200 EOFs retained).
-Parameter errors are constant over the current assimilation cycle.
-The parameter error prior pdf is Gaussian: -The prior pdf is kept unchanged between assimilation cycles.
Step c: augmented control vector prior probability distribution -Sample random parameter maps (100 members) from their pdf.
-Perform a model simulation for each member: the covariance of the ensemble forecast is the error covariance in the augmented control space.
-Parameters are constrained by bounds ⇒ the input parameter pdf is not really Gaussian.
-Assume a truncated Gaussian augmented vector distribution p(x) by imposing zero prior probability to large parameter increments.
Step d: parameter estimation using observations Apply observational update formulas to compute parameter corrections.
Atmospheric forcing parameters are the only source of error.
so that the prior pdf cannot be really Gaussian.In practice, this means that each time that maps of parameter increments are sampled from the prior distribution N (0, P α ) and added to the reference parameter maps, all values falling outside the physical bounds are reset to the value of the closest bound.
For instance, a negative precipitation value is reset to 0, or a cloud fraction value exceeding 1 is reset to 1.This set of operations implicitly define the prior pdf that is effectively assumed.This also means that the parameter perturbation may not be constant in time (and consequently that the correction may not be exactly piecewise constant) as soon as parameter values are found outside of their physical bounds.
As a last approximation in our experiments, the prior pdf for the error on the parameters is kept unchanged from one assimilation cycle to the next.This means that it is assumed that nothing is learnt about the parameters of the current cycle from the previous estimates.This is quite an important difference with respect to the experiments performed by Skachko et al. (2009), with the advantage that we do not need to parameterize the time dependence of parameter errors (since zero correlation is assumed).It is also safer to keep a zero mean error pdf around the reference parameter value.In that way, we can be certain to avoid any drift of the parameters from the reference (as observed in Skachko et al., 2009) and we do not need to add feedback to the reference (as they did) to prevent for filter instability.

Augmented control vector prior probability distribution
In order to use the ocean observations to estimate the parameters, we need to derive a prior probability distribution for the augmented control vector (step c in Table 1), including the ocean state and the forcing parameters (as explained in Sect.2.2).In our experiments, this prior probability distribution is approximated by a 100-member sample, that is obtained by sampling the forcing parameter probability distribution N (0, P α ) (described in Sect.3.3): α (i) , i=1, . . ., n=100 (using the method described in the appendix of Fukumori, 2002) and by performing the corresponding ensemble model forecast for the current assimilation cycle: x (i) , i=1, . . ., n=100.The 100 model forecasts, with their associated parameter maps ] represent the sample that we need to parameterize the prior probability distribution for the augmented control vector.This sample characterizes the sensitivity of the model forecast to parameter error around the current state of the system.This is a very important point because this sensitivity depends very much on the current ocean state.This is why it is impossible to perform the ensemble forecast once for all; it must be computed for each assimilation cycle from the current initial condition.Moreover, since the forcing parameters are the only source of error that is introduced in our ensemble experiments (no perturbation of the initial condition, no other model noise), the resulting probability distribution represents only that part of the total error that is caused by the forcing parameters.This is fully consistent with the experimental setup described in Sect.3.1, since we only seek to control the forcing parameters, but this also means that all other sources of errors in the system must be considered as observational error (and thus included in the parameterization of the observation error covariance matrix).
From the 100-member ensemble forecast x(i) , we parameterize the prior probability distribution of the augmented control vector as a Gaussian distribution N (x b , P), where xb is the background forecast obtained with zero parameter perturbation, and P is given by In parameterizing N (x b , P), we do not use the mean and covariance of the sample as xb and P because model nonlinearities can create a bias between x b and the sample mean x (i) .In our experiments, we want to improve the background model solution x b , not the ensemble mean x, which is never used as best estimate of the state of the system (as it could be in ensemble methods).To be consistent, we thus also need to characterize the sensitivity of the model forecast around x b and not around x.An additional reason is that the value obtained for the bias x−x b is related to the shape of the prior pdf for the forcing parameters p(α), which is not likely to be very accurately represented by our arbitrary Gaussian choice N (0, P α ).This bias problem arises because we try to solve a non-Gaussian problem approximately using a Gaussian approach.A rigourous solution can thus only be obtained by moving to a more general non-Gaussian scheme for the observational update.In this paper, we choose to leave these developments for further studies and to use the above Gaussian parameterization as an approximation.
Moreover, as an additional approximation, we only retain the first 50 principal components of the covariance matrix defined by Eq. ( 19), representing in general most of the total sample variance (around xb ). Figure 5 presents one column of the resulting correlation matrix (correlation with respect to SST at 66 • E, 1 • 52 S), for T a , q a , C and C E , showing for instance that the correlation is the largest close to the reference SST location, and decreases with the distance (as a general behaviour).It is dominantly positive for air temperature T a and negative for the other parameters q a , C and C E , consistently with the common physical sense.The horizontal shape www.ocean-sci.net/5/403/2009/Ocean Sci., 5, 403-419, 2009 of the correlation structure is highly anisotropic, with zonal correlations (along the equator) remaining significant over larger distance than meridional correlations, as a direct consequence of the anisotropy of the equatorial dynamics.We can even identify the correlated and anti-correlated separation zones for q a and C, to the separation between the North Indian Ocean currents (influenced by the Asiatic Monsoon) and the currents of the South Indian Ocean (influenced by the atmosphere anticyclonic circulation).The figure also shows that, due to the low rank (r=50) parameterization of the covariance matrix P, the correlations do not vanish at long distances as they do in the real world.This is why, in order to compensate for this deficiency in the parameterization of P, we impose vanishing long range correlation coefficients (see Testut et al., 2003;Brankart et al., 2003), using local observational updates (still performed in a reduced dimension space using Eq. ( 17), with locally defined Ŝ and R matrices).
After the observational update, it may be that the updated parameters do not satisfy the inequality constraints given by Eq. ( 18).If this situation occurs, the out-of-range parameter values are simply reset to the closest valid value as explained in Sect.3.3 for the ensemble experiments.This is of course an additional and quite crude approximation in the computa-tion of the posterior parameter estimates.The difficulty originates from the assumption of a constant parameter increment that is added to parameter maps that are not constant over the assimilation cycle.It is thus impossible to impose inequality constraints on the increment, and to use them to improve the shape of its prior probability distribution (for instance, by a truncated Gaussian assumption, as in Lauvernet et al., 2009 or by a nonlinear change of variable, as in Béal et al., 2009).

Truncation of the prior Gaussian distribution
In any inference problem, the accuracy of the posterior estimates crucially depends on the quality of the assumptions on the prior probability distributions, which in our problem are parameterized as Gaussian distributions: p(x)=N (x b , P), p(y|x)=N ( Ĥx, R).In the first distribution, we only include error that are due to forcing parameters, so that all other sources of errors must be included in the second distribution to correctly explain the dispersion of the observations (i.e.their covariance must be included in the observation error covariance matrix R).In our experiments, underestimating R is dangerous, because this means giving too much confidence to the observations, or in other words, interpreting an excessive part of the misfit with respect to the observations as due to forcing errors.The direct consequence is an excessive correction applied to the forcing parameters, that corresponds to very low prior parameters probability.Prohibitive values of the parameters, never occuring in the real system, can be reached because of the excessive tendency of fitting the observations (whose dispersion can only be explained by the existence of other sources of error in the system).Naturally, overestimating R is also dangerous, because this means not exploiting enough the observational information, and missing a part of the error variance that can be explained by forcing errors.
In order to reconcile the necessity to maintain the estimated forcing parameters inside a realistic range of values with the difficulty of producing a parameterization of the observation error covariance matrix R that is sufficiently accurate, we decide to proceed in the following way.First, we use a quite crude parameterization for the observation error covariance matrix: uncorrelated errors (diagonal R matrix) with uniform and quite small standard deviation: 0.1 • C for SST observations ans 0.02 psu for SSS observations.We can thus be quite sure that too much confidence is given to the observations (underestimated R).But second, we truncate the prior probability p(x) by imposing zero prior probability to large parameter increments.More precisely, the distribution in the reduced space p(ξ ) is truncated by the inequality constraints: |ξ i |≤γ , i=1, . . ., r. |ξ i | values larger than γ are thus assumed impossible.This corresponds to excluding an increment of the parameters along each error mode (each eigenmode in Eq. ( 14), left equation) if it is larger than γ times the standard deviation along that error mode.In our experiments, we set γ =3, which excludes any increment (along each mode) that is outside the 99.7% prior Gaussian confidence interval (i.e.occcuring typically in 0.3% of the param-eter maps sampled from a free model simulation, according to a Gaussian assumption).
From this modified prior probability distribution p(ξ ), it can be deduced from the Bayes theorem (8) that the posterior pdf p(ξ |y) is given by the same solution (17) as in the Gaussian problem but truncated by the constraints |ξ i | ≤ γ (since there is only multiplications by zero in Eq. ( 8), see Lauvernet et al., 2009, for more detail).The difference is that the previous Gaussian best estimate ξ a may not satisfy the constraints, and thus no more corresponds to maximum probability (but to zero probability).With the set of simple constraints |ξ i |≤γ , it is not difficult to see that the new maximum probability is obtained for Since the r inference problems are still independent, the maximum joint probability is indeed obtained if each onedimensional pdf is maximal (i.e.nearest to ξ a i within the valid interval).With this assumption, we can thus compute from ξ a i the forcing parameters that are capable of explaining the largest part of the misfit with respect to the observations, while remaining in a realistic range of variation (along each of the error modes).In that way, we expect that we can answer to the initial question (Sect.3.1): identifying what part of the misfit between model simulation and reanalysis can be explained by errors in the atmospheric forcing function.

1-year model simulation with parameter optimization
In this section the results of the free simulation are compared to those obtained with parameter optimization (without direct correction of the model state).In Fig. 1, the red line corresponds to the SST or SSS time evolution of the RMS error (difference with respect to the reanalysis) in the optimal simulation, that includes the application of our correction of the atmospheric forcing parameters.Starting from the perfect initial condition as indicated in Sect.3.2, we observe in the first assimilation cycles that the initial error trend is strongly reduced by the parameter correction (compare with the black curve starting from the same condition).This already shows the ability of the scheme to control a significant part of the model error that is due to the original forcing parameters.The large error reduction observed at the end of the initialization period then stabilizes over time over the one year diagnostic period.Much of the SST and SSS error in the free simulation is cancelled and the error variances (over the full domain) becomes about 3 times smaller than the corresponding values in the free simulation.Both curves stabilizes in time at around 0.5 • C RMS for SST and around 0.15 psu RMS for SSS with a drift that is now almost fully under control.However, this error remains very inhomogenous horizontally.Figure 6 shows maps of the spatial distribution of the SST and SSS misfit in terms of systematic error (top panels) and error standard deviation (bottom panels).As in Fig. 2, both statistics are computed for the one-year diagnostic period.Globally, by comparison to Fig. 2, we observe an important reduction of the SST and SSS systematic error and standard deviation everywhere in the global ocean.The only regions where a significant residual error remains are the Gulf Stream and Kuroshio regions, the Confluence region and the ACC and, to a lesser degree the Eastern Pacific equator (for SST) and the Western Pacific equator (for SSS).These errors are the consequence of the presence of other sources of error in the system which it is impossible to correct by just optimizing the forcing function.In particular, in the Western boundary currents and the ACC, an important part of the original error is due to the bad representation and localization of the ocean currents in our low resolution ocean model.This is why, in this case the parameter correction corresponding to an important bias is, for a large part, unrealistic as will be explained in Sect.4.3.
In order to have a better view of what occurs in the other regions (covering most of the ocean), the results in Fig. 1 are also presented (in the right panel) by excluding from the average the 10% of the ocean with the largest free run RMS misfit (for SST and SSS respectively).As compared to the full ocean results (left panel), the error reduction is here even more significant, with RMS misfits stabilizing at about 0.25 • C for SST ans 0.01 psu for SSS (without any residual drift).It appears that without considering the problematic areas listed above, the optimization of the atmospheric forcing parameters has a significant positive impact on the simula-tion, leading to surface ocean properties that are in very good agreement with the Mercator-Ocean reanalysis, and this result concerns up to 90% of the ocean surface.

Diagnostic of the mixed layer properties
In order to provide an idea of the vertical structure of the RMS error, on temperature and salinity, Fig. 7 shows the misfit with respect to the reanalysis for the simulation obtained with (dotted lines) and without (solid lines) parameters optimization.The figure is organized according to zonal bands: the Northern zone (between 55 • N and 19 • N, top panels), the Tropical zone (between 19 • N and 22 • S, middle panels) and the Southern zone (between 22 • S and 56 • S, bottom panels).Each panel shows the results for the Atlantic ocean (black), the Pacific ocean (red) and the Indian ocean (green).In general, dotted lines show smaller RMS misfit than solid lines, indicating the positive impact of the parameters optimization over the whole depth of the mixed layer.However, this impact is more significant in the mid-latitudes than in the tropics because the large subsurface difference observed at the equator, are associated to the thermocline (and thus well below the mixed layer), and are not connected to errors in the heat and fresh water flux.They correspond to errors in the depth of thermocline that result from wind forcing differences (likely connected to errors in the wind stress) or from the effect of data assimilation in the reanalysis.On the contrary, in the mid-latitudes, the error becomes much more constant along the vertical if the atmospheric parameters are optimized.The large error close to the surface in the free simulation (without parameter optimization) corresponds indeed to errors in the mixed layer that are clearly due to forcing errors, and that can be very substantially reduced by the optimization of the forcing parameters.

Diagnostic of the parameter estimates
In the previous sections, we have analyzed the impact of the parameter optimization on the temperature and salinity fields, and demonstrated globally that it produces thermohaline properties of the mixed layer that are in much better agreement with the Mercator-Ocean reanalysis.However, these positive results concerning temperature and salinity do not mean necessarily that the parameters themselves have been improved.This is much more difficult to demonstrate because our experiments are not twin experiments, so that we do not know the true values of the parameters.This is very different from the previous study by Skachko et al. (2009) who used twin experiments to demonstrate the accuracy of the parameter estimates.In our experiments, only indirect arguments can be proposed to study the relevance of the corrected atmospheric parameters.This can be done by trying to detect the two situations in which the T/S fields can be improved by irrelevant parameter corrections: -the optimization scheme produces irrelevant parameter corrections that compensate for other sources of error, -the optimization scheme produces irrelevant parameter corrections that compensate each other.
The first situation means that perhaps too much confidence has been given to the observations, and that part of the innovation is unduly attributed to atmospheric parameter errors.
In this, we must also include compensations for wind errors, which have not been included in the control vectors of the assimilation scheme (see Sect. 3.3).And the second situation means that the parameters may not be simultaneously controllable by the observations, i.e. the problem may be underdetermined.
The evaluation of the relevance of our parameter corrections will be based on two diagnostics: the time average of the parameter increment over the diagnostic period (Fig. 8) and the time standard deviation of the parameter increment (Fig. 9).The average increment can be compared to the mean parameter map in the original dataset (in Fig. 3) to get an idea of the average relative correction.And the standard deviation of the increment can be compared to the standard deviation of the time variability in the original data (in Fig. 4) to get an idea of the importance of the corrections with respect to the natural variability of the parameters.On these maps, we can see that almost identical corrections are computed everywhere for C E and C H consistently with their modelling by the aerodynamic bulk formulas.(Both are linear functions of the turbulent friction velocity.)Since the perturbations applied to the parameters in the ensemble forecast have the same covariance as a free model simulation, C E and C H can only be corrected in that way by the assimilation scheme.The first important thing to notice is that the amplitude of the correction is never larger than a few times the standard deviation of the natural variability of the parameters.This is the direct consequence of the limitation γ =3 that we have imposed on the correction in the reduced space.Moreover, approaching this maximal correction in average (i.e. a factor larger than 1 and even approaching 3 between Figs. 4  and 8) means that the correction is saturating on the γ =3 limit, which already indicates that the scheme is attempting to correct large temperature and salinity misfits that cannot be fully explained by the postulated level of parameter inaccuracy.This occurs first in regions of strong TS fronts resulting from important currents, where the scheme compensates the misplacement of the currents by irrelevant corrections of the parameters.In the Gulf Stream region, for instance, saturated corrections of the cloud fraction and relative humidity are produced North of the real position of the front (as represented by the reanalysis) to compensate for the overshooting of the current in the model solution.The same phenomenon happens in the ACC, but over a much larger area and with an even stronger impact on the parameters (see the mean and standard deviation of the cloud fraction and relative humidity around 60 • S).A similar problem also occurs in the Equatorial regions where surface temperature differences (certainly induced by wind errors) are here compensated by heat flux corrections.(The negative average precipitation increment in the Western Pacific is applied by the scheme to compensate for the negative salinity bias in this region; compare Figs. 2  and 6).In all these problematic regions, the saturation of the parameter increment (with the γ mechanism) also explains why it is also in these regions that SST and SSS differences with respect to the reanalysis (shown in Fig. 6) remain the largest.There, the scheme refused the large parameter corrections that would have been necessary to fully compensate the SST/SSS differences.The consequence is that the large SST/SSS misfits remain, while the parameters stay inside a reasonable range.
As a distinct kind of problem, it is interesting to remark the very large correction applied in the Southern ocean (along the Antarctic coast) to the C E and C H coefficients on the one hand, and on the air temperature T a on the other hand.These very large increments are not there to compensate very large SST and SSS errors (compare Figs. 2 and 6), so that they mainly compensate each other to produce the required SST and SSS corrections.This behaviour denotes the difficulty to control simultaneously several parameters using only SST and SSS observations.In this particular case, the large corrections are made possible by the very large standard deviation of these parameters in the natural variability of this region, and the scheme exploits this possibility as much as possible to fit the observations.There, the covariance of the parameter variability is certainly inadequate to represent correctly parameter errors.
The previous discussion summarizes the list of regions where the limitations of the method are most obvious.But everywhere else, the corrections are much smaller than the error bars that have been imposed and nevertheless sufficient to obtain temperature and salinity improvements (as observed by comparing Fig. 6 to Fig. 2).Such a result demonstrates that for the largest part of the ocean (in absence of other sources of errors strongly influencing SST and SSS), it is possible to produce moderate atmospheric parameter corrections that can drive the mixed layer properties of longterm ocean simulations very close to reanalysis data.

Influence of initial condition errors
Our first concern was to investigate ways of correcting the errors due to the atmospheric forcing and to dissociate them from other sources of error, like intrinsic model error or initial conditions errors.Up to here, we focused on that by starting the simulations from a perfect initial condition.As an additional experiment it is however useful to study the influence of initial condition error.For that purpose, we started a new free model simulation with parameters optimization from the initial date (30 December 1992) of the original reference simulation (without the relaxation that was performed to reach the perfect initial condition).This simulation (still without state correction) is illustrated by the green line in Fig. 1.What we observe first is the rapid error decrease during the first assimilation cycles of the experiment, showing the ability of the scheme to reduce the SST and SSS error that is present in the initial condition and to control the model error due to original parameters forcing.The comparison of the corresponding SST and SSS misfits with those obtained with perfect initial conditions (red lines), shows that the two experiments have the same kind of asymptotic behaviour on the long term, which means that the initial condition is progressively forgotten with time.Even if there is an impact of the initial error on the long term behaviour, particularly obvious for SSS, both simulations are characterized by error with similar magnitude over the diagnostic period, which means that the method can be applied with a similar success in presence of initial condition errors.

Conclusions
In this study, data assimilation experiments have been performed with the aim of controlling the parameters governing the atmospheric forcing, using idealized sea surface temperature and salinity observations, that are extracted from the Mercator-Ocean reanalysis data set.The results show that it is possible to compute piecewise constant parameter corrections, with predefined amplitude limitations, so that long-term free model simulations become much closer to the www.ocean-sci.net/5/403/2009/Ocean Sci., 5, 403-419, 2009 C. Skandrani et al.: Controlling atmospheric forcing parameters of global ocean models reanalysis data, with misfit variances typically divided by a factor 3 (for the global ocean) or by a factor 5 (if we exclude the frontal zones).However, the model that is used to perform these experiments is a low resolution model that does not represent correctly the Western boundary currents and other important circulation features which depend on resolution.The consequence is that part of this model error is incorrectly ascribed to the parameters so that the prescribed amplitude limitations saturate in these regions, thus indicating that the parameter corrections are irrealistic.Such problems can only be circumvented either by improving the model (for instance by increasing the resolution) or by controlling this error by data assimilation (for instance using altimetric observations).On the other hand, our experiments also suggest that a large part of the error (i.e. the misfit with respect to the reanalysis) can be explained by a bias on the reference parameters, with the consequence that our estimation scheme cannot be considered optimal (since centered prior probability distributions are assumed).All these results point towards the need for accurately specifying the prior parameter probability distribution and, despite of the deficiencies just mentioned, the experiments performed in this study already represent a significant step in this direction: by constructing the prior distributions locally in time, and by imposing strict limitations to the amplitude of the correction, we can be sure at least (by construction) that the parameter estimates always remain in a realistic range of values (i.e.inside their local range of variation in the input atmospheric data).From a methodological point of view, the application of such state dependent prior constraints is made practically possible by the use of a Monte Carlo method to simulate the joint parameter/state probability distribution.For that purpose, a truncated Gaussian assumption is used to parameterize these distributions, so that the posterior parameter estimates can be computed very efficiently.However, in the present study, the constraints have been defined according to a statistical criterion (99% confidence interval), which is certainly not the best way of summarizing the prior information about the parameter range of validity.In order to improve the definition of the estimation problem, the best perspective is certainly to give a more physical basis to the specification of the constraints.This can only be done by defining directly the range of validity in parameter space (and no more in the reduced space as in this paper), so that the simplifications that we brought to the truncated Gaussian filter of Lauvernet et al. (2009) would no more be applicable.Even if the algorithm can become somewhat more complex, this potential solution is certainly a good candidate to continue improving the prior parameter probability distribution, which we have just shown to be a key issue in the computation of more realistic parameter estimates.

FFig. 1 .Fig. 1 .
Fig. 1.Misfits between the simulations and the SST (top panels) or SSS (bottom panels) observations for the world ocean (south of 70 • N) as computed, without masking any regions (left panels) or by masking the 10% of the ocean surface that is charcterized by the largest free run RMS misfit (right panels).The figures show the free simulation (upper black solid line), the relaxation towards a perfect initial condition (dashed black line), the modified free simulation starting from this perfect initial condition (lower black solid line), the simulation with parameter optimization (green line), the simulation with parameter optimization starting from perfect initial condition (red line).The vertical dashed-dotted line marks the beginning of the one year diagnostic period.

Fig. 3 .
Fig. 3. Mean parameters in the free model simulation over the period 1992-1998.The figure shows C E ×10 −3 (top left panel), C H ×10 −3 (top right panel), C (middle left panel), P (in mm, middle right panel), T a (in K, bottom left panel) and q a (bottom right panel).

Fig. 4 .
Fig. 4. Parameters standard deviation in the free model simulation over the period 1992-1998.The figure shows C E ×10 −3 (top left panel), C H ×10 −3 (top right panel), C (middle left panel), P (in mm, middle right panel), T a (in K, bottom left panel) and q a (bottom right panel).

Fig. 5 .
Fig. 5. Correlation with respect to SST at 66 • E, 1 • 52 S, for T a (top left panel), q a (top right panel), C (bottom left panel), C E (bottom right panel).

Fig. 6 .
Fig. 6.Maps of systematic error (top panels) and error of standard deviation (bottom panels) for SST (in • C, left panels) and SSS (in psu, right panels) as obtained for the model simulation with parameter optimization (starting from perfect initial condition).

Fig. 7 . 30 Fig. 7 .
Fig. 7. RMS misfit with respect to the reanalysis temperature (left panel) and salinity (right panel).The Figure is organized according to zonal regions: the Northern zone (between 55 • N and 19 • N, top), the Tropical zone (between 19 • N and 22 • S, middle) and the Southern zone (between 22 • S and 56 • S, bottom).Each panel shows the results for the Atlantic ocean (black), the Pacific ocean (red) and the Indian ocean (green).Solid lines correspond to the free simulation and dashed lines correspond to the simulation with optimized atmospheric parameters.For this figure, the 10% of the ocean with the largest free run RMS misfit have also been excluded from the computation of the horizontal averages.

Fig. 8 .
Fig. 8. Time average of the optimized parameter increment over the diagnostic period.The figure shows this result for parameters C E ×10 −3 (top left panel), C H ×10 −3 (top right panel), C (middle left panel), P (in mm, middle right panel), T a (in K, bottom left panel) and q a (bottom right panel).

Fig. 9 .
Fig. 9. Standard deviation of the optimized parameter increment over the diagnostic period.The figure shows this result for the parameters C E ×10 −3 (top left panel), C H ×10 −3 (top right panel), C (middle left panel), P (in mm, middle right panel), T a (in • K, bottom left panel) and q a (bottom right panel).