Turbulent viscosity optimized by data assimilation

As an alternative approach to classical turbulence modelling using a first or second order closure, the data assimilation method of optimal control is applied to estimate a time and space-dependent turbulent viscosity in a three-dimensional oceanic circulation model. The optimal control method, described for a 3-D primitive equation model, involves the minimization of a cost function that quantifies the discrepancies between the simulations and the observations. An iterative algorithm is obtained via the adjoint model resolution. In a first experiment, a k ± L model is used to simulate the one-dimensional development of inertial oscillations resulting from a wind stress at the sea surface and with the presence of a halocline. These results are used as synthetic observations to be assimilated. The turbulent viscosity is then recovered without the k + L closure, even with sparse and noisy observations. The problems of controllability and of the dimensions of the control are then discussed. A second experiment consists of a two-dimensional schematic simulation. A 2-D turbulent viscosity field is estimated from data on the initial and final states of a coastal upwelling event.


Introduction
In ocean models, the system of mean equations for momentum, temperature and salinity contain unknown second order correlations between¯uctuations components of momentum and the buoyancy. Closure assumptions need to be introduced to relate frictional eects to the calculated large-scale velocity ®eld, and the diusive eects to the gradients of the temperature and salinity ®elds. These eects are produced in turbulence-like fashion by motions on scales too small to be resolved by the grid. These are called sub-grid scale phenomena. To close the system, the eects of these scales are usually parametrized in a simple way through the use of the eddy viscosity and diusivity concept. Typically in the past, constant values of eddy viscosity (with dierent values in the horizontal and vertical directions) or empirical variable eddy viscosity, allowing non-linear eects to be more realistic, were used. More recently, with the increasing resolution capacities of present-day computers, the diusion coecients have often been determined from local values of scalar properties of the turbulence. Approaches of this type make use of one or more additional transport equations (Rodi, 1980). Today, while models have become quite sophisticated in computational and turbulence modelling aspects, they depend on numerous calibration parameters that must be extrapolated from limited ®eld measurements. Furthermore, many uncertainties that cannot be accounted for by these turbulence models remain.
As an alternative approach, an inverse strategy is here proposed. In applied sciences, inverse methods are commonly used to extract useful inferences about the world from physical measurements (Menke, 1984). Here, the turbulent viscosity and diusivity are ®tted from data issuing from oceanographic observations. The method is a speci®c case of the optimal control method often called the adjoint method (Seiler, 1993). The ideas and mathematical concepts of optimal control theory were formalized about thirty years ago (Lions, 1971) and have since received much attention for applications in oceanography (Begis and Crepon, 1975;Devenon, 1990). The aim of the optimal control method is to ®nd the best parameters of a model to simulate the computed values closest to those observed. The adjoint method has been most often used to ®t initial conditions (Moore, 1991), but also boundary conditions , or speci®c parameters like the ocean surface heat uxes (Roquet et al., 1993)). This variational method involves the minimization of a cost function which is the norm of the dierence between the computed and observed values. An algorithm is obtained via the adjoint equations for the construction of the gradient of the cost function with respect to the parameters. Once the gradient has been determined, the minimization can be performed using any gradient descent algorithm.
The eciency of the optimal control method to optimize time-constant viscosity distributions in 1-D vertical models has already been underlined (Yu and O'Brien, 1991;Panchang and Richardson, 1993). In these studies, initial and boundary conditions are assumed to be well known, the variational formulation of the dynamic model acting as a strong constraint. In the case of 1-D advection-diusion equation (Leredde et al., 1998) it has been shown that, with such a strong constraint formulation, the simultaneous optimization of the initial conditions, the boundary conditions and the space and time distributed viscosity coecient becomes rather dicult. This is due to the fact that a single solution for all the controls is not ensured if no penalization of the ®rst guesses of the controls is included in the cost function. Instead of using this type of strong constraint formulation, Eknes and Evensen (1997) have developed a weak constraint formulation which allows the observational data, the turbulent vertical viscosity, the initial and boundary conditions and also the model equations to be aected by errors. This approach involves a sophisticated iterative representer method (Bennett, 1992) and generates a highly non-linear problem whose solution depends on the required ®rst guesses of the control and the error weights. Eknes and Evensen (1997) have used a 1-D Ekman model which is only an approximation of the primitive equations. Their results indicate that model de®ciencies, such as neglected physics, are accounted for through the weak constraint formulation to ensure an inverse solution in agreement with the measurements. In this respect and in the framework of this non-exact model, their ®t is closer to the data than those obtained by the strong constraint formulation (Yu and O'Brien, 1991).
In the present study, a 3-D primitive equation model has been adopted. Since the model is expected to be more exact and the representer method is not yet implemented in this three-dimensional case, the strong constraint variational formulation is used. Furthermore, in previous studies (Yu and O'Brien, 1991;Panchang and Richardson, 1993;Eknes and Evensen, 1997), the turbulent viscosity was taken to be constant in time. In fact, as mentioned many times in the literature (e.g. Blumberg and Goodrich, 1990), turbulent mixing is time-dependent especially for transient phenomena such as wind-driven or tidal processes. Therefore, the optimal control method is formalized in the more general frame of the optimization of a turbulent viscosity which can vary in the four space and time dimensions.
This induces a large number of parameters to be optimized requiring at least the same number of data to be assimilated. We have thus chosen to create synthetic observations with a complete turbulence model. The chosen model is a k v model (Leendertse et al., 1973;Leendertse and Liu, 1977) using a single transport equation for the turbulent kinetic energy k and an algebraic formulation for a mixing length scale v. It allows the construction of a data set and of the associated turbulent viscosity set to be recovered by the inverse method.
The physical and mathematical background is given in Sect. 2, including the description of the 3-D multilayer dynamic model and the adjoint model resulting from the variational formalism. This section deals with the complete and continuous optimization problem and adopts a more general framework than the applications which are made in Sects. 3 and 4. In addition to the numerical computational cost, it becomes rather dicult to solve optimization problems which have a size equal to the number of grid points and time steps in the discretized form. Therefore, the ®rst applications of our data assimilation model concern fairly simple oceanographic processes allowing the numerical treatment of the method with an acceptable space-time dimension.
In Sect. 3, a wind-driven mixing event of a strati®ed water column inducing inertial oscillations is simulated by the 3-D multilayer dynamical model with open lateral boundaries and consequently the process is 1-D. Section 4 deals with the optimization of a 2-D space-dependent turbulent viscosity in the case of a wind-driven mixing event of coastal upwelling. Concluding remarks are presented in Sect. 5.

The circulation model
The direct model (Leendertse et al., 1973), recently implemented (Lellouche, 1995) is a classical primitive equation model for coastal seas. Many others are based on a similar approach (e.g. Nihoul, 1977;Thouvenin and Salomon, 1984;Van Dam and Louwersheimer, 1990). The Boussinesq approximation (density is constant except in the buoyancy term, (Pedlosky, 1982)), and those of hydrostasy and f -plane are made. The physical state of the¯uid is given by the seven equations for the average variables: In this set of equations, uY vY w are the mean velocity components in the three space dimensions xY yY zY t denotes the time, p the pressure, the salinity, the temperature, q the sea-water density. q 0 Y 0 Y 0 are constant reference values for density, temperature and salinity. f is the Coriolis parameter, g the gravitational acceleration, a the saline contraction coecient, a the thermal expansion coecient. m h Y m t are the horizontal and vertical turbulent viscosity coecients. m h Y m t Y m h and m t are the turbulent diusivity coecients. A ®rst approach, the zero order turbulence closure, consist of giving constant values to these viscosity and diusivity coecients without using additional transport equations. When needed, heat¯uxes at the surface and bottom and surface stresses can be introduced. In the following, only wind stress and bottom stress are considered and modelled.
At the upper boundary, the horizontal wind stress vectors can be estimated by the quadratic expression: where q is the air density, w the wind speed vector at 10 m above mean sea surface and g d a drag coecient. The bottom stress vectors is also expressed in the same way: where q is the water density at the bottom layer, is the horizontal velocity vector at the bottom layer and g the Chezy coecient, homogeneous with g 1 2 . For the numerical implementation (Lellouche, 1995), a ®nite dierence approximation has been used both in time and space according to a vertical multilevel geometrical discretization.
In a ®rst model version, used for open boundary condition optimizations , the turbulent viscosity and diusivity coecients were all set to constant values. A simple study of sensitivity of the model to the horizontal viscosity and diusivity coecients values shows that the horizontal diusion terms are helpful to damp down the small-scale instabilities. Nevertheless a ®rst approximation given by constant horizontal coecients seems to be sucient. The horizontal diusion terms in Eqs. (1), (2), (5) and (6) were already written under this hypothesis.
In contrast, the vertical turbulent viscosity and diusivities have a real in¯uence on the numerical results, particularly for strati®ed seas, and must be speci®ed with more accuracy.

The k + L turbulence model
A classical approach consists in using a turbulence model where the turbulent viscosity and diusivity coecients are expressed as functions of the turbulent state of the sea water. This turbulent viscosity and diusivity concept gives rise to various models. Among these models, the length scale formulation (Blackadar, 1962) can be chosen and might be corrected to take correctly into account the strati®cation eects. This dependence can be introduced via the Richardson number (Munk and Anderson, 1948). To go further, second order formulations (Rodi, 1980;Mellor and Yamada, 1982), known as k À e and k À v models, have became rather popular, notably in the coastal and estuaries oceanographic community, and satisfy operational numerical needs. This operational goal restricts the implementation of higher order formulations with transport equations for Reynolds stresses (Mellor and Yamada, 1974;Andre et al., 1979). In coastal zones, recent works (Nihoul et al., 1989. Luyten et al., 1996 have shown that cruder models with only one transport equation for the turbulent kinetic energy (k-models) could be suciently ecient. In the framework of this study and in order to progress towards our assimilation aims, a k-model, initially developed by Leendertse and Liu (1977), has been chosen. This model, called k v model, takes the following form.
The vertical eddy viscosity coecient m t is evaluated by using a buoyancy extended Prandtl-Kolmogorov hypothesis: where v is a mixing length scale, k the turbulent kinetic energy, and m and g m positive constants. i is the``turbulent'' (not the``gradient'') Richardson number.
The mixing length scale is given by an algebraic function of the distance to the bottom where j is the Von KaÂ rmaÂ n constant and d the local water column height. The vertical eddy diusivities for heat and salinity transport are assumed to be proportional to the turbulent vertical viscosity: where rY are respectively the turbulent Prandtl and Schmidt number. The turbulent kinetic energy is governed by the following transport equation: where m k h , the horizontal turbulent diusivity for k, is considered as constant, and m k t , the vertical turbulent diusivity for k, is assumed to be proportional to the turbulent viscosity: where r k is also often called the turbulent Prandtl number for the turbulent kinetic energy. Inside the¯uid, the shear production of turbulent energy is expressed by: The bottom layer friction is taken into account by an additional production term: where g is the Chezy coecient previously introduced.
In the top layer, the wind and waves eects give also rise to an additional term: where a is a non dimensionless m À1 positive constant. The buoyant production q is modelled according to Instead of solving a transport equation which involves more model assumptions than the corresponding equation for kY e, the turbulent kinetic energy dissipation is modelled according to where e 0 is a positive constant.

The optimal control formulation
Even if the applications presented in Sects. 3 and 4 are two dimensional (1D space and 1D time or 2D space) and in order to adopt a more general framework, the complete 4D formalism is here described. The optimal control aims at ®nding the best parameters of the model to simulate computed values closest to those observed. It is exempli®ed here when the turbulent viscosity coecient set is considered as the control to be optimized. The physical state is described by the vector uY vY wY pY Y Y q, solution of the set of dierential Eqs (1)±(7) which can be re-written: where p is a non-linear operator, X is the space domain and 0Y s is the time interval. X is considered in this general development as three dimensional. Here, m m t Y m t Y m t is the only control of the system. In a more general case, the control could be the initial conditions (Le Dimet and Talagrand, 1986), the boundary conditions (Lellouche et al., 1994) or other tunable parameters such as friction coecients (Begis and Crepon, 1975), optimized separately or together (Devenon, 1990;Leredde et al., 1998). In the following, in order to simplify somewhat the problem the turbulent Prandtl and Schmidt numbers, respectively r and , are set to constant values. Without restraining the scope of the present study, taking r 1, the turbulent diusivities for heat and salt transport will be equal to the turbulent vertical viscosity. The control is then restricted to m t as a space and time dependent variable. In this general framework, m t is a function of the three space dimension xY yY z and the time dimension t.
The observations of being designated as , a cost function which measures the mean squared discrepancy between the model solution and the observations can be de®ned as: where`Y b X Â 0Ys is the inner product in v 2 X Â 0Y s, this space being de®ned as the set of square integrable functions over X Â 0Y sX g is the covariance operator associated with the observations. Penalty terms could be added to this functional formulation. The convergence of the following method is then improved by the recall to a``®rst guess'' of the control (Yu and O'Brien, 1991) or a smoothing term of the distributed control (Panchang and Richardson, 1993).
The optimal control m t opt minimizes the cost function t . The optimal control set cannot be obtained directly. It must be reached by successive minimization of t through a descent algorithm. A classical way consists in computing the gradient of t relative to the control. Introducing the adjoint state Ã (Lions, 1971), the expression of the gradient is: and the adjoint state veri®es that: where dj denotes the adjoint operator. The details of the method for deducing the adjoint model can be found in Leredde et al. (1998) for the case of a one dimensional Burger's equation and in Leredde (1999) for these primitive equations and will not be presented further here. From Eq. (24), using the proper de®nition for p r 2 y i , is the variance of the observation y i corresponding to the variable y i in the case of statistical independence between data measurements. In fact, r 2 y i , is the inverse of the weigh given to each observation. For example, if no temperature data is available, the weight is zero, r 2 is taken to be a large number and so there is no forcing term in Eq. (30).
Keeping this notation y i for any of the seven state variables, the dimensionless cost function can be written.
The expression of the gradient is then given by: If the turbulent viscosity is not fully time and space dependent, this general expression is integrated on the domain where the turbulent viscosity is constant. In Sect. 3, m t is considered as horizontally constant and the temperature is not taken into account. The expression of the gradient of the cost function relative to m t zY t becomes where C is the horizontal domain. In Sect. 4, m t is a function of two space dimensions x and z. The only assimilated data are temperature ones and the expression of the gradient of the cost function relative to m t xY z is where v is the domain of variation of y. In fact, the general expression (33) can be adapted for each new con®guration. For example, without the salinity term and with a time integration, the gradient expression of previous studies (Yu and O'Brien, 1991;Panchang and Richardson, 1993;Eknes and Evensen, 1997) could be found.
Once the gradient has been computed, the minimization of the cost function can be performed using any gradient descent algorithm. Here, a classical quasi-Newton method (Gillbert and Lemarechal, 1989) is used.
For the numerical implementation, this continuous formulation must be converted to the discretized one. The way to proceed is not straightforward and involves theoretical and technical developments .

Description of the experiment
Although the numerical tools described in the previous section should allow the treatment of a fully three dimensional case of a rotating¯uid, a simple physical 1-D case is investigated as a ®rst attempt. The 3-D model is then used with open lateral boundaries at which Neumann conditions are imposed and all the variables are taken to be horizontally uniform. In this academic framework, the observation data set is created using the k v model described in Sect. 2.2. This numerical construction of error free data is usually done to check the results of new assimilation methods (Courtier and Talagrand, 1990). The aim is to recover these``simulated'' data by means of the full identi®cation of the turbulent viscosity coecients. It means that a zero order turbulence closure version of the model is used instead of solving the k v closure.
In order to illustrate this 1-D optimization problem, the turbulent mixing of the water column resulting from a wind energy input at the sea surface is studied. The direct simulation consists in modelling the development of inertial oscillations in the presence of a halocline break. An 8-day simulation is performed with a time step of 20 min and a vertical grid size of 2 m. The model is started from a state of rest considering a simpli®ed two-layer water column with fresh water (10 m thick, 0 psu) upon salty water (20 m thick, 30 psu). The wind blows along the x direction and its speed increases linearly from 5 m s À1 to 10 m s À1 during the simulation. The solution of the computation is given in Fig. 1. Inertial oscillations remain trapped in the upper layer of the water column until the mechanical mixing energy, resulting from the downward diusion of a fraction of the energy input by wind, overcomes the stabilizing eect of the strati®cation. Then, the halocline is suddenly broken and the motion can diuse downward and the salinity upward. In the following, these simulated data are used as observations to be assimilated by the optimal control method. An initial model solution (Fig. 2) from which to start the minimizing procedure is computed with a constant value of 10 cm 2 s À1 arbitrarily chosen for the turbulent viscosity. One can remark that this value is not accounted for as a ®rst guess in the cost function as it is needed for a weak constraint formulation (Eknes and Evensen, 1997) to ensure a single solution (Bennett and Miller, 1991). In a strong constraint formulation, even though it can be used to accelerate the convergence of the method (Yu and O'Brien, 1991), it is not required if the problem remains over-determined (Menke, 1984). In any case, such a ®rst guess choice could be dicult in this case of timedependent viscosity.
A comparison between the observation simulation ( Fig. 1) and the initial simulation (Fig. 2) shows that the solution of the model is highly controlled by the turbulent viscosity. This latter is about zero at the halocline that acts as a physical wall for the diusion. As the wind speed increases, the eddy viscosity increases to a threshold value inducing the sudden halocline break. The oscillation amplitude and damping rate follow this time-dependence of the viscosity. As this viscosity varies in the whole space and time domain X Â 0Y s, the dimension of the control vector is equal to w Â x where w is the number of grid points and x is the number of time steps. In order to keep the problem over-determined, the number of data to be assimilated must overtake the control dimension. All the model solution variables, here uY vY can be taken as observations. Dierent numerical experiments using dierent data sets are performed. The variances of the selected data are then set to a ®nite value corresponding to the inverse weights given to each observation. The other data variances are set to in®nity expressing the complete lack of any other available information. The forcing Fig. 1a±d. Solution of the direct simulation used to construct the observations. a Velocity ucm Á s À1 , b velocity vcm Á s À1 , c salinity (psu), d turbulent viscosity cm 2 Á s À1 . The contour intervals are 5 units for each variable plotted terms which appear in the second members of Eqs. (25)± (31) are then all equal to zero except for those corresponding to the chosen observations.

Optimization results
The most basic numerical experiment is to consider all the assimilated data uY vY as error free data. In comparison with a cost function formulation including a ®rst guess (Yu and O'Brien, 1991), the decrease of the cost function t is slower, but after about 3000 iterations, t reaches almost zero. The optimized viscosity values (Fig. 3) provide ®elds of uY vY which appear to be indistinguishable from the observations. This results shows the feasibility of the method without using a ®rst guess. Nevertheless, it must be recalled that such a good ®t between model solution and data set is possible only because the data set was obtained via the complete k v model and so is fully consistent with the optimization constraints. This cannot be expected with in-situ experimental data.
Besides, the optimal set of turbulent viscosity values (Fig. 3) is sometimes dierent locally from the one given by the k v model used to construct the data set. It can be noticed that, even when dierent, the turbulent viscosity issued from the physical turbulence model and the optimal viscosity issued from the data assimilation ®tting can lead to the same numerical solution of the model. This paradoxical behaviour can be explained. The optimal viscosity has a degree of freedom (DoF) per grid point and time set while m t , computed with the k v model, has a smaller number of DoF equal to the number of constants appearing in the k v turbulence model. In the optimization procedure, the turbulent viscosity value at a grid point and time step is not Fig. 2a±c. Solution of the initial simulation with a constant turbulent viscosity m t 10 cm 2 Á s À1 . a Velocity ucm Á s À1 , b velocity vcm Á s À1 , c salinity (psu). The contour intervals are 5 units for each variable plotted Fig. 3. Optimal turbulent viscosity cm 2 Á s À1 . The contour intervals are 5 cm 2 Á s À1 constrained by the global physical consistency due to the equation of the k v turbulence model. Furthermore, the sensitivity of a model to its parametrization may be so weak that the control of these parameters cannot be ensured. As already explained by Panchang and Richardson (1993), the eddy viscosity appear only in the local diusion terms so that in weak vertical gradient regions, it can take any value without in¯uencing the simulation results. This states the problem of the possibility to control models by parameters that do not have everywhere a real in¯uence on the solution.
The ability of the method to recover the model control from a part of the observation data set must now be tested. In a second experiment, where only one variable, for instance , is assimilated, the optimization procedure succeeded in recovering exactly the salinity ®eld t 0 and the velocity ®elds (Fig. 4) with a better accuracy in strati®ed regions. The controllability problem leads to an optimal eddy viscosity that can take any possible value in zones of zero vertical salinity gradient. The uniqueness of the solution is then not ensured and the turbulent viscosity remains near its initial value here 10 cm 2 s À1 . For practical purposes, the velocity ®eld is then aected in these regions and it could be worse with another initial value like 0 cm 2 s À1 or 100 cm 2 s À1 . The retrieval of the velocity ®eld from the salinity ®eld is interesting in itself since hydrological data are often more easily available than dynamic data. It must however be borne in mind that in the model the salt diusivity is assumed to be proportional to the turbulent viscosity for the momentum diusion ( is a constant). Further generalization to models with eddy viscosity and diusivity using other relationships should be made.
For the same reason, the assimilation of only one velocity component u or v leads to similar conclusions. u and v are in quadrature, and so give the same content of information for the optimization. The sudden event of the halocline break is included in the velocity ®eld information. In fact, with a slower convergence rate compared with the previous experiment, the model solution is recovered even for the salinity.
The three previous numerical experiments were performed with a great number of DoF with good eciency due to the ideal consistency of the data with the model. Unfortunately, oceanic observations are often sparse and noisy. In order to have a more realistic representation of in-situ observations, the following assimilation experiments are carried out with noisy and sparse data. To simulate dierent station measurements conducted from an oceanographic vessel, either hydrological or current pro®ling experi- Fig. 4a±d. Optimal solution and optimal turbulent viscosity obtained by the assimilation of the salinity data. a Velocity ucm Á s À1 , b velocity vcm Á s À1 , c salinity (psu), d turbulent viscosity cm 2 Á s À1 . The contour intervals are 5 units for each variable plotted ments, the proposed application consists in assimilating noisy velocity and salinity pro®les taken at 48 h, 96 h, 144 h and 192 h. A white and gaussian noise is added with a rms equal to 6 cm s À1 for the velocities and 4 psu for the salinity. This rather severe noise has been chosen to test the robustness of the optimal control method. As the control dimension greatly oversizes the number of sparse data, the unique determination of an optimal time-dependent eddy viscosity cannot be expected. To keep the problem over-determined, the control dimension must be reduced at least to the same order as the number of data. Following the idea of Panchang and Richardson (1993), a penalty term e.g. dm t tadt could be added in the cost function to keep a continuous time-dependence of the viscosity. Here, a time-constant turbulent viscosity pro®le over each 48-h time-interval is chosen. The problem consists then in optimizing four viscosity pro®les given for each period: 0 h±48 h, 48 h±96 h, 96 h±144 h, and 144 h±192 h. The pro®le assimilation leads to the recovery of the full solution uY vY zY t (Fig. 5). One can note the time jumps of the solution corresponding to the optimal turbulent viscosity pro®les. Figures 6 shows the uY vY optimal pro®les at t 96 h (the other time step results are not shown). Even if the assimilated pro®les are noisy, the optimal pro®les are very close to the true pro®les issued from the full solution given by k v model. The initial starting pro®les are computed with a constant viscosity value 10 cm 2 s À1 . Even with severe noise, the model is able to extract the mean information and to supply results provided with some physical consistency given by the model equations. The optimal eddy viscosity pro®les (Fig. 7) evolve with time and with the wind increase. It is almost zero on the halocline for the two ®rst pro®les (48 h, 96 h) before it breaks. When the turbulent viscosity does not control the solution, its value remains near its initial starting value. The non-sequential nature of the proposed assimilation scheme must be underlined. Indeed, the information over the entire data recording is used as a whole to optimize the turbulent viscosity value at a grid point and time step. A sequential optimization would conversely consist in the separate processing of juxtaposed time periods with initial and ®nal state observations. In the actual method, observations on past and future states are accounted for in addition to the initial and ®nal states of the considered sequence. The case of a wrong data pro®le at a given time, say at 48 h, which exhibits for instance a poor vertical resolution leading to an anticipated halocline break, could be corrected with later pro®les where this discontinuity remains.
To summarise the previously described numerical experiment, it can be said that the assimilation technique turns out to oer a satisfactory robustness with noisy and sparse synthetic data. Fig. 5a±c. Optimal solution obtained by the assimilation of sequential and noised data (48 h, 96 h, 144 h and 192 h). a Velocity ucm Á s À1 , b velocity vcm Á s À1 , c salinity (psu). The contour intervals are 5 units for each variable plotted

A two space dimension optimization case
The aim of this experiment is to show that the optimal control method can be used for more than one space dimension. The method described in Sect. 2.3 in the more general case has been applied in the previous section for only one space dimension and the time dimension. The control size was already very large in this previous simple case, involving a high numerical cost and a slow convergence of the method. This cost is enhanced by the need of the direct solution to solve the adjoint model: the values of the direct variables are to be stored in the direct path and read from the memory in the inverse path. This is why the control is here limited only to two space dimensions and considered as constant in time. In fact, a space dimension takes the place of the time dimension. For the same reason, a simpli®ed theoretical model of ocean and of wind patterns are considered using a simple geometrical grid. The integration proceeds also for a relatively short simulation time. In addition it is the temperature distribution that here modi®es the strati®cation.
The 2-D space-dependent viscosity optimization problem is exempli®ed on the simulation of a transient coastal upwelling assuming that m t is time-independent. The domain is closed to the west by a coast on which zero cross-shore transport conditions are speci®ed. A constant wind 5 m s À1 is assumed to blow towards the north. By assuming the domain is open to the south and to the north, the solution does not vary along the direction y. By taking Neumann open ocean boundary conditions to the south and north, a two space dimension problem in a plane perpendicular to the coast is obtained. Neumann conditions are also imposed at the east boundary. The domain is extended to 40 km from the coast with a horizontal grid size Dx 2 km, whereas the depth is constant (40 m) with a vertical grid size Dz 2 m. The latitude is 45 N, which leads to a Coriolis parameter f 10 À4 s À1 . The time step is Dt 50 s for a 48 h simulation. The model is started from a state of rest, a surface temperature of 25 C and a linear vertical gradient of temperature D Dz 0X5 C m À1 . Salinity is taken to be uniform and constant in time.
The observations are created from a 2-D steady viscosity ®eld corresponding to the results given by the k v model at t 48 h. At this time, the turbulent viscosity ®eld is nearly steady and so can be taken as the steady control to be recovered. Nevertheless, one can remark that, even if the viscosity is considered as constant in time, coastal upwelling phenomena is already not steady after 48 h. During this transitional period, inertial waves always propagate (see Crepon andRichez, 1982, or Kundu et l., 1983, for more details). The model therefore gives results as expected (e.g. Foo, 1981). Close to the coast, nearly 5 km, a coastal jet expands in the direction of the wind. The surface mixing layer is advected towards the open sea, whereas the return current stays between the boundaries of the bottom mixing layer. Near the coast, cold waters are advected from the bottom to the surface. Figure 8 presents the solution uY vY wY Y m t at the end of the observations simulation.
The aim is then to reconstruct the 2D space-dependent turbulent viscosity ®eld from the solution of the model. The initial conditions (state of rest and linear temperature strati®cation) are presumed to be known. The initial solution obtained with an uniform viscosity value of 10 cm 2 s À1 (Fig. 9) is considered as the starting point of the optimization procedure. The results obtained for this initial simulation show that it is necessary to consider a spatial variability of the viscosity.
In order to elaborate the problem and make the method more adapted to usual practical oceanographic application, only temperature data are assimilated, assuming the velocity ®eld to be unknown. As the Fig. 8a±d. Solution of the direct simulation used to construct the observation of coastal upwelling at t 48 hX x (km) is the distance to the coast. a Velocity vcm Á s À1 , b Temperature C), c vertical current vector uY w Â 10 3 , the maximum length is 19 cm s À1 , d turbulent viscosity cm 2 Á s À1 control is considered as constant in time, a single time sequence of data is needed.
First, it has been chosen to assimilate the complete and noise-free temperature ®eld at t 48 h. The number of data is almost equal to the number of parameters and it is already a signi®cant diculty to suppose that no information is known for the velocity ®eld and that no other data are available between 0 h and 48 h.
The optimal control method provides by inversion a viscosity ®eld which leads to a quite accurate simulation of the temperature ®eld. One can remark that, as well as with the salinity data in the previous example, the only assimilation of the temperature data allows to reconstruct the dynamics. Figure 10 presents this complete solution, very close to the simulated observations. Secondly, in order to test the robustness of the method, the previous experiment is realised assimilating the complete but noise-distributed temperature ®eld at t 48 h. A white and gaussian noise is added to the temperature ®eld. The level of noise is raised step by step as long as the method remains ecient enough. The deduced borderline case, presented here, corresponds to a noise with a rms equal to 2 C, which is more higher than the in-situ measurement errors occurring with usual devices, such as CTD pro®lers. Figure 11 shows the assimilated temperature ®eld.
It can be seen that the added noise is smoothed by the physics of the model which acts as a strong constraint. The discrepancies of the data with the model, induced by a gaussian and white noise, remains weak enough to allow the method to work as shown in Fig. 12. This 2-D optimization case will allow us to consider more complex situations in the future.

Summary and conclusions
As an alternative to classical turbulence modelling, the data assimilation method of optimal control is formulated to optimize the turbulent viscosity in a 3-D primitive equation model of the ocean circulation. The chosen model and its full turbulence formulation k v is described and used to simulate synthetic observations (velocity, salinity or temperature) of two wind-driven events: the ®rst concerns an open strati®ed ocean, and the second, a coastal upwelling. These data are assimilated to reconstruct the turbulent viscosity ®eld, which depends either on time and 1-D space or 2-D space. An Fig. 9a±c. Solution of the initial simulation of a coastal upwelling with a constant turbulent viscosity m t 10 cm 2 Á s À1 at t 48 hX x (km) is the distance to the coast. a velocity vcm Á s À1 , b temperature ( C), c vertical current vector uY w Â 10 3 , the maximum length is 17 cm s À1 adjoint model of the direct model (without the complete turbulence formulation) is integrated backward in time to compute the gradient of the cost function with respect to the control and therefore performs the minimization.
The interest of ®tting the vertical turbulent mixing via m t must be stressed. In many coupled physicobiogeochemical models (e.g. Pinazo et l., 1996), the turbulent state input is only described by a distribution of the eddy vertical viscosity. Its optimization by data assimilation could be useful in this case, the detailed understanding of the turbulent processes not being required. Of course, the success of the method is conditioned by the nature of the data. In the current work, it is pointed out that a vertical mixing process can be recovered using only salinity or temperature data. Furthermore, it would also be interesting to apply this method to reconstruct a physical mixing event from the observations of a conservative tracer (Vandenberghe, 1992).
The optimal control method has shown its eciency in reconstructing the turbulent viscosity ®eld from synthetic velocity or salinity data. Next, as has already been done on a time-independent 1-D optimization in an Ekman model (Yu and O'Brien, 1991), the optimization in a 3-D primitive equation model will have to be tested with assimilation of real observations. This preliminary study has pointed out the diculties that could be encountered. Fig. 10a±d. Optimal solution and optimal turbulent viscosity obtained by the assimilation of the temperature data of the coastal upwelling. x (km) is the distance to the coast. a Velocity vcm Á s À1 , b Temperature ( C), c vertical current vector uY w Â 10 3 , the maximum length is 14 cm s À1 , d turbulent viscosity cm 2 Á s À1 Fig. 11. Assimilated and noisy temperature data (°C) of the coastal upwelling, at t 48 h, obtained by addition of a white and gaussian noise with an rms equal to 2 C The ®rst problem is the large dimension of the control vector. This involves an important numerical cost and a slow convergence rate of the method. Then, a number of data oversizing the control dimension is required. The eddy viscosity variations should be chosen versus the available data. For example, in Sect. 4, as only initial and ®nal states are known, the eddy viscosity is optimized as a time-constant. The full time-space 4-D variational formalism has been described and should be adapted for each physical situation and assimilable data.
The second problem is the controllability problem. Prior knowledge of the distributed turbulent viscosity is sometimes not available, especially for a time and space dependent problem, and the method must be investigated without ®rst guess. As the eddy viscosity is not always an in¯uential parameter in the calculations, its uniqueness may not be ensured in some areas. Uniqueness in these zones could be achieved using a ®rst guess but with a loss of accuracy of the model results in the strong eddy viscosity in¯uenced regions. Moreover, the optimal viscosity can then take any values which appear sometimes to be physically inconsistent. In the optimization procedure the zero order turbulence closure model allows a degree of freedom by grid point and time step. In contrast, in the k v turbulence model, the viscosity is strongly constrained by the k and v via Eq. (10), and so has a smaller number of degrees of freedom equal to the number of constants e.g. g m Y kY aY e 0 . In this case, turbulence model calibration is required. It could be viewed as a challenging problem to consider in turn these parameters as controls of the model to be optimized. However, hardly any diculties are probably to be expected due to the strong non linearities of the dependence between the model solution and these constants. Fig. 12a±d. Optimal solution and optimal turbulent viscosity obtained by the assimilation of the noisy temperature data of the coastal upwelling. x (km) is the distance to the coast. a Velocity vcm Á s À1 , b temperature ( C), c vertical current vector uY w Â 10 3 , the maximum length is 14 cm s À1 , d turbulent viscosity cm 2 Á s À1