Estimation of friction parameters in gravity currents by data assimilation in a model hierarchy

Abstract. This paper is the last in a series of three investigating the friction laws and their parametrisation in idealised gravity currents in a rotating frame. Results on the dynamics of a gravity current (Wirth, 2009) and on the estimation of friction laws by data assimilation (Wirth and Verron, 2008) are combined to estimate the friction parameters and discriminate between friction laws in non-hydrostatic numerical simulations of gravity current dynamics, using data assimilation and a reduced gravity shallow water model. I demonstrate, that friction parameters and laws in gravity currents can be estimated using data assimilation. The results clearly show that friction follows a linear Rayleigh law for small Reynolds numbers and the estimated value agrees well with the analytical value obtained for non-accelerating Ekman layers. A significant and sudden departure towards a quadratic drag law at an Ekman layer based Reynolds number of around 800 is shown, in agreement with classical laboratory experiments. The drag coefficient obtained compares well to friction values over smooth surfaces. I show that data assimilation can be used to determine friction parameters and discriminate between friction laws and that it is a powerful tool in systematically connecting models within a model hierarchy.


Introduction
The realism todays and tomorrows numerical models of the ocean dynamics is and will be governed by the accuracy of the parametrisations of the processes not explicitly resolved and resolvable in these models.A typical example is the thermohaline circulation of the world ocean, a basin-scale Correspondence to: A. Wirth (achim.wirth@hmg.inpg.fr)long-time circulation.It is key to the climate dynamics of our planet and it is governed by small scale processes.The thermohaline circulation starts with the convective descent of dense water masses at high latitudes, a process performed by plumes smaller than a kilometre (Wirth andBarnier, 2006, 2008).When the dense water touches the ocean floor and/or passes through a strait, a gravity current is created, which is strongly influenced by the friction forces at the ocean floor (Wirth, 2009).The bottom friction is influenced by the dynamics at scales of the order of one meter and less.Similar small scale processes are likely to determine the subsequent deep western boundary current and the upward transport and the mixing in the interior and coastal regions of the ocean.These small scale processes will not be explicitly resolved in ocean general circulation models (OGCMS) even in a far future.The understanding and parametrisation of these processes is of paramount importance to the progress in modelling the dynamics of the climate.
In the present work I focus on bottom friction determined by small scale three dimensional turbulence.This process has much smaller spacial scales and faster time scales than the large scale circulation above.This scale separation in space and time is the prerequisite for an efficient parametrisation of bottom friction.Two types of bottom friction laws are commonly employed in engineering and geophysical fluid dynamics applications: a Rayleigh friction with a friction force per unit area F = ρhτ u depending linearly on the fluid speed u at some distance of the boundary and a quadratic drag law F = ρc D |u|u, where τ and c D are the linear friction parameter and the drag coefficient, respectively, h is the thickness of the fluid layer and ρ the density of the fluid.In todays ocean general circulation models a mixture of both laws is commonly employed, where the friction force per unit area is given by (1) Published by Copernicus Publications on behalf of the European Geosciences Union.
A. Wirth: Estimation of friction parameters in gravity currents by data assimilation The velocity c represents unresolved velocity due to tidal motion and other unresolved short time-scale processes.Typical values for c are a few tenths of centimetres per second.Such friction law leads asymptotically to a linear friction law for small velocities u c and to a quadratic drag law for u c.
The precise determination of friction laws and parameters are also fundamental to the understanding of the large scale dynamics in the ocean.Grianik et al. (2004), showed that the energy-containing scales of the large scale ocean dynamics is determined by the friction laws and parameters.
Bottom friction has also been identified as an important process acting as a sink of kinetic energy.Kinetic energy is principally injected into the ocean by the surface wind-stress at the basin scale.It can only be dissipated at scales where molecular viscosity is acting, that is below the centimetre scale.The pathway of the energy from the basin scale to the dissipation scale is currently to large parts unexplored.The turbulent bottom boundary layer dissipates energy at a rate proportional to c D V 3 where c D is the (local) drag coefficient and V the (local) flow speed near the boundary.The fact that the energy dissipation is proportional to the third power of the speed emphasises the importance of high speed events and processes in the vicinity of the ocean floor.The determination of the precise value of the drag coefficient which varies over an order of magnitude depending on the roughness of the boundary (see e.g., Stull, 1988) is, therefore, key to determining the energy fluxes and budget of the worlds ocean.
In the present work I focus on the dynamics of oceanic gravity currents which is governed by bottom friction, as it was shown in Wirth (2009).Gravity currents are clearly high-speed events at the ocean floor with speeds of over 1 ms −1 .Friction laws and parametrisations at solid boundaries have been a major focus of research in fluid dynamics for the last century, due to their paramount importance in engineering applications.Although large progress has been made in determining the friction over rough boundaries (see e.g., Jiménez, 2004) the friction depends of a variety of properties of the ocean floor which are undetermined, as for example: the roughness type, the roughness scales, the multiscale properties of its roughness, the sediment suspension, the orientation of the roughness elements and the variability of the roughness, to mention only a few.Furthermore, these properties will not be available to ocean modellers in the foreseeable future.In the present work I determine the friction laws and parameters by observing the time evolution of the thickness of a gravity current.In the absence of a sufficient density of data from oceanic observations or laboratory experiments, the data is taken from non-hydrostatic numerical simulations of gravity current dynamics.Once the friction laws and parameters determined, they can be used in future integrations of ocean general circulation models (OGCM).To solve the inverse problem of determining the friction laws and parameters from the evolution of the gravity current, I use data assimilation.The formalism of estimating friction laws in oceanic gravity currents has been introduced and discussed in Wirth and Verron (2008).The purpose of the present work is not simply to use data assimilation so that the shallow water model mimics the non hydrostatic dynamics, but to extract the friction parameters and laws from the non hydrostatic dynamics.
I am here interested in the dynamics that governs the thick part of the gravity current, called the vein (see Wirth, 2009).The friction layer, the thin part at the down-slope side of the vein, represents a water mass that is lost for the gravity current, does not contribute to its further evolution, and is mixed into the surrounding water in a short time.The dynamics of the friction layer, which is less than 20 m thick, is also likely to be determined by small scale structures of the ocean floor.
This paper is the last in a series of three.In the first (Wirth and Verron, 2008) we determined the feasibility and convergence of estimating friction parameters while observing only the thickness of the gravity current.In the second (Wirth, 2009), I studied the dynamics of an idealised gravity current using a non-hydrostatic numerical model.I showed that such a gravity current has a two part structure consisting of a vein, the thick part, which is close to a geostrophic dynamics perturbed principally by the influence of the Ekman pumping caused by bottom friction.At the down-slope side of the vein is the friction layer, the other part, completely governed by frictional dynamics.In the present work the two approaches are combined.The data is taken from non-hydrostatic numerical simulations of Wirth (2009) and provided to the assimilation scheme introduced in Wirth and Verron (2008).The reader interested in the detailed description of the data assimilation is invited to consult the latter and the reader eager to know more about the dynamics of the gravity current and its turbulent characteristics should consult the former paper.For the present paper to be self-contained, an overlap to the two previous papers is unavoidable.
The methodology presented here, although developed for the case of idealised gravity currents, is not restricted to such, but can be generalised to other processes investigated by a model hierarchy.Conceptually, I show how to systematically connect models of different complexity in a model hierarchy.The same methodology can be employed when parameters are to be estimated using data from observations or laboratory experiments.An example is given in a series of papers by Ghil (1997), Sun et al. (2002) and Kondrashov et al. (2008), where parameters in an idealised coupled oceanatmosphere model are estimated using an extended Kalman filter.This methodology was further applied to a highly nonlinear problem in solid mechanics by Kao et al. (2006).Parameter estimation using the extended Kalman filter in an intermediate complexity earth system model are given in Hargreaves et al. (2004) and Annan et al. (2005).They adjust the climatology of the ocean and atmosphere model to observed data by estimating O(10) parameters and subsequently perform simulations of climate change scenarios.In the next section I introduce the physical problem considered and the two mathematical models employed to study its dynamics, followed by a discussion of their numerical implementations.The data assimilation algorithm connecting the two models is discussed in Sect.3. A detailed presentation of the experiments performed is given in Sect.4, results are presented in Sect. 5 and discussed in Sect.6.
2 Idealised oceanic gravity current on the f-plane

The physical problem considered
In the numerical experiments I use an idealised geometry, considering an infinite gravity current in a rotating frame on an inclined plane with a constant slope, and I do not allow for variations in the long-stream direction.As discussed in Wirth (2009) a gravity current dominated by rotation is, to leading order, in a geostrophic equilibrium where the downslope acceleration due to gravity is balanced by the Coriolis force.Such gravity current flows along-slope, not changing its depth (see Fig. 1).It is friction that makes a rotating gravity current flow downslope.This is the opposite in nonrotating gravity currents, where friction opposes the downslope movement and results from studies of non-rotating gravity currents can not be applied to rotating gravity currents.In the geometry considered here, I neglect the long-stream variation of the gravity current.Such a dynamics is usually referred to as 2.5-D as it includes the fully three dimensional velocity vector but depends on only two space dimensions.Please note, that such simplified geometry inhibits large scale instability, the formation of the large cyclones and other large-scale features, which is beneficial to our goal of studying the friction laws due to only small scale dynamics.
The initial condition is a temperature anomaly which has a parabolic shape which is 200 m high and 20 km large at the bottom.The velocity is initially geostrophically adjusted (see Sect. 2.2 for details).
The physical problem is considered with the help of two mathematical models of different complexity.The first are the Navier-Stokes equations with a no-slip boundary condition on the ocean floor, subject to the Boussinesq approximation.The second is a single layer reduced gravity shallow water model.The bottom friction is explicitly resolved in the first while it has to be parametrised in the second.

The Navier-Stokes model
The mathematical model for the gravity current dynamics are the Navier-Stokes equations subject to the Boussinesq approximation in a rotating frame with a buoyant scalar (temperature).The state vector is formed by the temperature anomaly of the gravity current water with respect to the surrounding water and all the three components of the velocity vector.I neglect variations in the stream wise (y-) direction of all the variables.Such type of model is usually referred to as 2.5-D.The x-direction is upslope (see Fig. 1).This leads to a four dimensional state vector depending on two space and the time variable, ( T (x,z,t),u(x,z,t),v(x,z,t),w(x,z,t)).
The domain is an inclined rectangular box that spans 51.2 km in the x-direction and is 492 m deep (z-direction).On the bottom there is a no-slip and on the top a free-slip boundary condition.The horizontal boundary conditions are periodic.The thickness of the gravity current is denoted by the variable h(x,t).The initial condition h(x,t = 0) is a temperature anomaly which has a parabolic shape which is 200 m high and 20 km large at the bottom, as described in the previous subsection.The magnitude of the temperature anomaly T is varied between the experiments.The reduced gravity is • T with g = 9.8066 ms −2 .The Coriolis parameter is f = 1.03 × 10 −4 s −1 .The initial velocities in the gravity current are geostrophically adjusted: The fluid outside the gravity current is initially at rest.The buoyancy force is represented by an acceleration of strength g cos(α) in the z-direction and g sin(α) in the negative xdirection to represent the slope of an angle α = 1 • .This geometry represents a rectangular box that is tilted by an angle of one degree.Such implementation of a sloping bottom simplifies the numerical implementation and allows for using powerful numerical methods (see next subsection).

Numerical implementation of the Navier-Stokes model
The numerical model used is HAROMOD (Wirth, 2004).HAROMOD is a pseudo spectral code, based on Fourier series in all the spatial dimensions, that solves the Navier-Stokes equations subject to the Boussinesq approximation, a no-slip boundary condition on the floor and a free-slip www.ocean-sci.net/7/245/2011/Ocean Sci., 7, 245-255, 2011 boundary condition at the rigid surface.The time stepping is a third-order low-storage Runge-Kutta scheme.A major difficulty in the numerical solution is due to the large anisotropy in the dynamics and the domain, which is roughly 100 times larger than deep.There are 896 points in the vertical direction.For a density anomaly larger than 0.75 K the horizontal resolution had to be increased from 512 to 2048 points (see Table 1), to avoid a pile up of small scale energy caused by an insufficient viscous dissipation range, leading to a thermalized dynamics at small scales as explained by Frisch et al. (2008).The horizontal viscosity is The anisotropy in the turbulent mixing coefficients reflects the strong anisotropy of the numerical grid.I checked that the results presented here show only a slight dependence on ν h ,κ h and κ v by doubling these constants in a control run.This is no surprise as the corresponding diffusion and friction times are larger than the integration time of the experiments.There is a strong dependence on ν v as it determines the thickness of the Ekman layer and the Ekman transport, which governs the dynamics of the gravity current as I have shown in Wirth (2009).The vertical extension of the Ekman layer is a few meters while the horizontal extension of the gravity current is up to fifty kilometres.In a fully turbulent gravity current the turbulent structures within the well mixed gravity current will be isotropic and will therefore measure only a few meters in size.To simulate a fully turbulent gravity current 10 5 grid points would be necessary in the horizontal direction to obtain an isotropic grid.This is far beyond our actual computer resources.
The time of integration is limited due to the spreading of the gravity current, the upper boundary of the gravity current is almost stationary the lower boundary slides down the incline.The rate of spreading depends on the initial density anomaly.The actual integration times in all experiments are given in Table 1.

The shallow water model
The second, and less involved, mathematical model for the gravity current dynamics is a 1.5-D reduced gravity (1.5 layer) shallow water model on an inclined plane.The shallow water model, first proposed by de Saint Venant (1871), and its various versions adapted for specific applications is one of the most widely used models in environmental and industrial fluid dynamics.For the derivation of the reduced gravity shallow water equations in a geophysical context, I refer the reader to the text book by Vallis (2006), and references therein.
As stated in the introduction I here specialise to a gravity current with no variation in the y-direction, the horizontal direction perpendicular to the down-slope direction.That is I have three scalar fields ũ(x,t), ṽ(x,t) and h(x,t) as a function of the two scalars x and t.The thickness of the gravity current is given by h.The velocities ũ(x,t)= h 0 u(x,z,t)dz/ h and ṽ(x,t)= h 0 v(x,z,t)dz/ h represent the vertical averages over the whole layer thickness h(x,t) of the local velocity components u(x,z,t) and v(x,z,t).I further suppose that the dynamics within the gravity current is well described by a two-layer structure: an Ekman layer with a thickness scale δ= √ 2ν v /f at the bottom and the rest of the gravity current above.To take into account several features of the gravity current dynamics discovered with the non-hydrostatic model (see Wirth, 2009), the shallow-water model employed in Wirth and Verron (2008) had to be refined.Results presented in Wirth (2009) show, in agreement with Ekman layer theory, that the x-component of the velocity is concentrated in the Ekman layer.This confinement of u to the Ekman layer is included in the shallow water equations through the Coriolis-Boussinesq variable: If there is no vertical shear β=1.If there is no bottom friction, the shear in a homogeneous fluid layer is small and β≈1.When Ekman layers are present there is substantial shear and β>1.Ekman layers are a conspicuous feature in all geophysical flows subject to bottom friction and this has to reflected in the mathematical model employed.In a shallow water model this can be done by explicitly via the Coriolis-Boussinesq variable.The Coriolis-Boussinesq variable is a function of the variation of the velocity profile in the vertical.When using the linear Ekman spiral, I obtain to first order β=h/(2δ) when h δ which depends on the layer thickness h and is thus a function of space and time.Please see the Appendix A for details concerning the calculations of the Coriolis-Boussinesq variable.The Coriolis-Boussinesq variable is widely used in mathematical modelling of river dynamics and appears naturally when the shallow-water equations are systematically derived for flow with vertical shear of the horizontal velocities (see Gerbeau and Perthame, 2001).
The governing equations for the velocities vertically averaged over the whole thickness of the dynamic layer h, Ekman layer plus the rest above, are given by: The left hand side terms include the reduced gravity g =g ρ/ρ the slope α and the Coriolis parameter f .The non-linear term in the ṽ-equation, that is ũ∂ x ṽ, is absent as the u velocity is concentrated in the Ekman layer and On the right hand side I have the terms involving dissipative processes.Please note that this includes the bottom friction as well as the friction at the interface, in numerical experiments with the non-hydrostatic model (Wirth, 2009) the latter is found to be smaller than the former.The parametrised friction is represented in the first term on the left hand side, involving: The Coriolis-Boussinesq variable β in Eq. ( 8) is explained in the Appendix.There are three free parameters τ,r and c D .There is a linear friction constant parametrising dissipative effects that can be represented by linear Rayleigh friction.This linear friction is represented by two parameters τ and r, the first represents the part that is independent of the thickness and the second is divided by the square of the thickness.Usually the thickness of laminar boundary layers grows in time, when the dynamics is influenced by rotation, this growth is halted creating the well known Ekman layer dynamics (see e.g., Vallis, 2006).Using only the term containing the parameter τ represents well the dynamics when an Ekman layer is developed, that is, when the thickness of the gravity current is larger than a few times the Ekman layer thickness.For smaller thicknesses, smaller than a few times the Ekman-layer thickness, a term involving the thickness of the gravity current should be included.I found heuristically that adding a friction force that depends quadratically on the layer thickness fits the data reasonably well for the gravity currents with a low density anomaly.In this publication I am, however mostly interested in the friction forces acting on the vein as stated in the introduction.In the results presented below the friction force due to the depth dependent term is a negligible part of the total friction, in the vein.
The quadratic friction drag c D , models the turbulent friction between the ground and the gravity current.The Coriolis-Boussinesq variable is given by β=h/(2δ) for large values of the layer thickness, as derived in Appendix A. If β is smaller than 2 it is put equal to 2 in the advection term and if it is smaller than 1/4 it is put to 1/4 in the friction terms.The first choice is consistent with the fact that in the advection term the Coriolis-Boussinesq variable is by definition always larger than unity.Concerning the friction terms: for large layer thickness, β agrees with the value calculated in the Appendix A and for a very small layer thickness, β=1/4 agrees with the fact that the friction in the x-and y-direction in Eqs. ( 8) and ( 9) should have the same form in the absence of an Ekman spiral.I checked that the choices for these thresholds have only a negligible influence on the results.
The term involving the viscosity/diffusivity ν represents horizontal dissipative processes, its value is chosen to provide numerical stability of the calculations (see Sect. 2.5).
Clearly, the geostrophic velocity 2: is a solution of Eqs. ( 4)-( 6) when D u =D v =0 and ν=0.The estimation of D or more precisely of the parameters τ , r and c D are the subject of the present work.By using data assimilation I plan to obtain the friction constants τ and c D and thus determine if the friction acting on the vein is dominated by a linear or a quadratic law.

Numerical implementation of the shallow water model
The shallow water model is implemented with a first order finite difference scheme in space and time.There are 500 points in the x-direction, leading to a resolution of 200 m, the time step is 5 s.The value of the horizontal viscosity/diffusivity is a function of the resolution and provides for the numerical stability of the code.I verified that the here presented results show only a negligible dependence when the value of the horizontal viscosity/diffusivity was halved and doubled, the actual value used is ν H =5 m 2 s −1 , it is identical to the corresponding value ν h in the non-hydrostatic model.

Ensemble Kalman filter and its implementation
For this publication to be self contained, the implementation of the Ensemble Kalman filter is explained here, the reader familiar with Wirth and Verron (2008) is invited to skip the section.Further details about the data assimilation procedure and the convergence properties are explained in (Wirth and Verron, 2008).The ensemble Kalman filter (EnKF) is the main tool of our experiments performing the parameter estimation and providing us with the actual parameter values.The EnKF was introduced by Evensen (1994) Brusdal et al., 2003).I refer the reader not familiar with the EnKF and the employed notation to the above mentioned publications.Every hour in time and every 1 km in the x-direction, the vertical extension of the gravity current, that is h(x,t), is assimilated.Choosing a horizontal resolution for the assimilation 5 times sparser than the dynamical model does not only reduce the size of the assimilation experiment, but is also consistent with the fact, that the grid-scale dynamics of the numerical model is dominated by dissipation and has only negligible dynamical information.I am only assimilating the vertical extension of the gravity current as it is the variable most easily measured in the ocean and in laboratory experiments.The measurement of the vertically integrated velocity within the gravity current, the other dynamical variable of the shallow water model, is more difficult to measure in the ocean and in laboratory experiments.
The assimilation is performed on the augmented state vector consisting of the vertical extension, the two velocity components and the three constant-in-time friction parameters: x( x, t) = h( x, t), ũ( x, t), ṽ( x, t),τ,r,c D t . (10) Where x and t is the discretized version of x at assimilation grid-points and t at assimilation times, respectively and t denotes transposition.The only (observed) variable assimilated is the vertical extension (h) of the gravity current.The analysis step for the Kalman filter (see Ghil et al., 1981 for an early geophysical application) reads, where the index i = 1,...,m runs over the realisations.The observation operator H projects the state vector in the space of observations.The noise vectors i represents the independent Gaussian-distributed zero-mean and σ -standarddeviation noise added to every observed value (see Burgers et al., 1998) and W=σ 2 I where I is the unity matrix.In the Extended Kalman Filter (EnKF) the estimated error covariance matrix is given by where m is the size of the ensemble and .denotes an ensemble average.The covariance matrix is truncated to a tridiagonal form to avoid spurious correlations between distant correlations, caused by under sampling.This is also consistent with the dynamics at hand as the maximum velocities (see Sect. 2.1) and wave speeds, given by g h are of the order of 0.2 ms −1 .This means that in the assimilation period of 1 h information can travel to the next assimilation point one kilometre away, but it can not reach over the distance of two assimilation grid points.In between assimilation points h is interpolated linearly.In general the observed value of the vertical extension of the gravity current h obs (x,t) includes measurement errors η(x,t) and is related to the true value by h obs (x,t)=h true (x,t)+η(x,t).For consistency the measurement error η(x,t) has the same first (zero mean) and second order moment ( σ 2 ) as the noise vectors i (x,t), but does not depend on the actual realisation, that is i.When assimilating data, σ has to be provided prior to the experiment whereas σ is usually not known and can only be estimated.
In our parameter estimation experiments, the ensemble size is m=100, this is much larger than the number of parameters to estimate, that is three, equal to the number of observations at each assimilation time, but smaller than the dimension of the augmented state vector x=(h, ũ, ṽ,τ,r,c D ), that is 303.Using an ensemble size an order of magnitude larger did not improve the convergence significantly, reducing the ensemble size an order of magnitude leads to a frequent divergence of the assimilation.
It is important to note that not only the parameter values are estimated, but the entire state vector is updated every time data is assimilated.This allows to perform parameter estimation in the case where the shallow-water model is unable, even with perfectly adjusted parameters, to reproduce some aspects of the Navier-Stokes dynamics.Entrainment is a typical example, it happens in the Navier-Stokes dynamics, there is no reliable parametrisation available for the shallowwater model and friction parameters can not account for all the impacts of entrainment on the evolution of the layerthickness.Trying to estimate the friction parameters without correcting the other values of the state vector frequently leads to a divergence of the estimation procedure in experiments where mixing, entrainment or detrainment is present (experiments with high values of the density anomaly).
Another important point is, that the parameter values are supposed to be constant in time.The values of the friction parameters are clearly non-negative, so every time the assimilation scheme provides a negative value of one of these parameters, which is possible due to the linearity of the analysis step and the statistical nature of the EnKF, the value is put to zero.

Experiments
The goal of our data assimilation experiment is to determine the friction laws acting on the vein of a gravity current, linear versus quadratic, and the corresponding friction parameters, by only observing the vertical extension (h) of the gravity current.The vertical extension, that is the density structure of a gravity current is the variable that is easiest to measure and to observe in the ocean and in laboratory experiments.Velocities, even their average values are hard to determine as they are highly variable in space and intermittent in time.
The thickness h is extracted from the runs of the Navier-Stokes model, it is determined by the position above ground of the thermocline that corresponds to half the initial maximum temperature anomaly in the gravity current at the assimilation points x, every 1 km, and the assimilation times t, every hour.It is the time series of these shapes h(x,t), and only these, which are provided to the assimilation run to determine the friction parameters.Please note that the space integral of h, the total volume, is not conserved but has some slight variation due to mixing, entrainment and detrainment, whereas in the shallow water model without assimilation it is conserved.The data assimilation procedure corrects the total state vector, including the h-value and the method can thus deal consistently with variations in the total gravity current fluid volume (or area in our 2.5-D case).
The important results concerning the data assimilation experiments, as for example: convergence properties of the friction parameters and a good choice of the observation frequency and other parameters are discussed in (Wirth and Verron, 2008).
A series of seven experiments was performed.The density anomaly T is the only parameter varied between the experiments.A list of all experiment performed, the duration of the experiment, the average geostrophic velocity vg and the Reynolds number, based on the Ekman layer thickness, given by δ= √ 2ν/f , are shown in Table 1.The Reynolds number based on the Ekman layer thickness, is the determining parameter for laminar flow.When vertical velocities develop (absent in the homogeneous linear Ekman layer), the now turbulent Ekman layer increases in thickness and the actual Reynolds number is larger than calculated by Eq. ( 14) and given in Table 1.Results shown in the next section show that vertical velocities become significant for temperature anomalies of about one Kelvin.When vertical velocities are important the flow is characterised by the surface Rossby number, where z 0 is the roughness length, which for smooth surfaces is given by z 0 =ν v /u * .The friction velocity u * = √ c D v g and the surface Rossby number are results of our numerical calculations, rather than an initial parameter, and is therefore presented in Table 2 of the results section.In the turbulent boundary layer vertical velocities of order u * develop.The typical length scale is now formed by a balance between rotation effects and vertical velocity δ turb ≈u * /f , rather than the balance between viscosity and rotation for the laminar case.Observations show that the actual thickness of the turbulent Ekman layer is ≈δ turb /4, including a log-layer of thickness ≈δ turb /10.The above teaches us, that the surface Rossby number is the equivalent to the Ekman-scale Reynolds number with a turbulent Ekman layer thickness δ turb replacing the laminar Ekman layer thickness δ.The former characterises turbulent flow and the latter laminar flow.
To avoid confusion I continue to discriminate between the experiments in terms of temperature anomaly rather than Reynolds or surface Rossby number.The surface Rossby number is also a measure of the thickness of the turbulent Ekman layer (or the log-layer) to the thickness of the viscous sub-layer δ sub =ν/u * (please see McWilliams, 2006, chapter 6 for a concise introduction to planetary boundary layer dynamics).I refer the reader to Wirth (2009) for a more detailed description of the gravity current dynamics considered here.
In the data assimilation experiment the observation error for the layer thickness was set to σ =10 m.I checked that the results presented in the next section do not depend significantly on the actual value of σ when chosen within reasonable limits.The assimilation is done every hour, the length of the experiment is given in Table 1 and assimilation interval in space is 1 km.Initially the width of the gravity current is 20 km, the experiments are stopped before it reaches 50 km.The ensemble of the initial values for the parameters was chosen randomly, with an ensemble mean value at the www.ocean-sci.net/7/245/2011/Ocean Sci., 7, 245-255, 2011 order of magnitude of the expected values.The procedure was iterated by using the ensemble at the end of a assimilation experiment as the initial ensemble for the consecutive experiment.I also tried inflation of the ensemble values, the distance from the ensemble mean was increased by a constant factor for every ensemble member at every assimilation time, to avoid spurious decrease of the ensemble variance.
Inflation did not change significantly the results.All results reported here were done without inflation.For further information concerning the assimilation scheme, its performance and convergence properties I like to refer the reader to Wirth and Verron (2008).
The non-hydrostatic calculations were done on 16 parallel processors at IDRIS high performance computation centre, whereas the assimilation of the data by the shallow water model where performed on a lap-top computer (please see Wirth, 2009;Wirth and Verron, 2008 for further details).

Results
In the case of a (linear stationary) Ekman dynamics, the linear Rayleigh friction can be calculated analytically, as already explained in Wirth (2009): with an Ekman layer thickness δ= √ 2ν/f .The value for the calculations presented here is: τ = 2.27×10 −4 ms −1 .
I started by performing parameter estimation experiments using the shallow water model used in Wirth and Verron (2009), that is, without the Coriolis-Boussinesq variable and without the thickness dependent parameter r, results were poor.The results for the vein improved significantly when the non-linear term was multiplied by the Coriolis-Boussinesq variable.A free evolution of the shallow water model (without data assimilation) with the estimated values of τ and c D showed results that compare poorly with the data for the thickness h from the non-hydrostatic calculations in friction layer, where the thickness is small.This problem improved when the parameter space was augmented by the thickness dependent parameter r.In the words of data assimilation: without the calculated Coriolis-Boussinesq variable and the estimated friction parameter r the assimilation converged to a (local) minimum of the cost function which has a value that is not small enough.In the words of a dynamicist: without the Coriolis-Boussinesq variable and the parameter r an important physical process is missing and the dynamics can not be understood without it.I like to emphasise that without using data assimilation one might have thought that one just did not find the right minimum in parameter space, with data assimilation one can be more confident, so not certain, that the minimum was found but the model had major deficiencies so that it is not capable to represent well enough the dynamics.Anyhow, adding the Coriolis-Boussinesq variable and the thickness dependent friction leads to an improvement of the representation of the dynamics, especially for the cases with a high temperature anomaly ( T ).
The experiments are listed in Table 1 and results are presented in Table 2.The results for the Rayleigh friction parameter (τ ) compare well to the analytical value of Eq. ( 16) for small values of the temperature anomaly.The impact of the thickness dependent friction parameter r is negligible in the thick part of the gravity current in all the experiments performed (see Table 2).The drag coefficient is small in these cases with a small temperature anomaly but increases abruptly with the temperature anomaly and represents a substantial part of the total friction in the experiments with a larger temperature anomaly.I also calculated for all experiments an effective drag coefficient: and the surface Rossby number Ro.It is the effective drag coefficient that is usually represented in diagrams like Fig. 2 (dashed line) by the engineering community (Schlichting and Gertsen, 2000).
The sudden appearance of a drag law at Reynolds numbers Re Ek ≈8×10 2 is shown.The values of the effective drag coefficient are shown in Fig. 2. Another interesting feature of the assimilation procedure is, that I not only obtain effective drag coefficient as is usually the case but I clearly manage to separate the friction process in a linear and a quadratic part, gaining further insight in the dynamics.
When comparing the dynamics of a free run of the shallow-water model, without data assimilation to the dynamics of the Navier-Stokes model, the agreement is generally satisfactory, but deteriorates with increasing temperature anomaly of the gravity current.I have to emphasise that a difference between a free run with adjusted parameters to the data, does not a priori devaluate the parameter estimation.The difference might originate from a process not represented, neither directly nor parametrised, in the shallowwater model.During the data assimilation this difference is corrected by the update of the whole state vector, in the free run this correction is absent.More precisely: entrainment, detrainment and mixing changes the layer thickness h and the overall volume of the gravity current, the importance of these processes increase with the Reynolds number.The corresponding parameters are not estimated in the shallow-water model.A change of the total volume can not be obtained by changing the friction parameters, as they do not affect the volume of the gravity current.The layer thickness is corrected at every data assimilation as it is part of the state vector.The free run and the data do therefore not necessarily agree, even when the total volume of the gravity current is considered.
As stated in Sect. 3 I supposed the parameters to be constant, a variability in time is clearly detected in all three parameters.Their variability is however small compared to their absolute value.Parameter values obtained by assimilating only the first half of the time series differed by less than 20% to assimilations with the total length.

Discussion
I have demonstrated, that friction parameters and laws in gravity currents can be estimated using data assimilation.The results clearly show that friction follows a linear Rayleigh law for small Reynolds numbers and the estimated friction coefficient agrees well with the analytical value obtained for non-accelerating Ekman layers.A significant and sudden departure towards a quadratic drag law at Reynolds number of around 800 is shown, roughly in agreement with laboratory experiments (Nikuradse, 1933;Schlichting and Gertsen, 2000).It is the effective drag coefficient that is usually represented in diagrams like Fig. 2 (dashed line) by the engineering community (Schlichting and Gertsen, 2000).When friction processes are considered experimentally, numerically or analytically the friction laws and coefficients depend on the Reynolds or surface Rossby number.The laws and coefficients are than compared in the different experiments and theories.Their overall behaviour consists always in a linear friction law for laminar flow (small Reynolds numbers) and a drag law (large Reynolds numbers) with a roughly constant drag coefficient for turbulent boundary layers.For intermediate Reynolds numbers it is difficult to separate the two friction laws without using data assimilation.Although the quadratic drag law is dominated by the linear friction at the Reynolds numbers considered in this study, the assimilation procedure clearly manages to detect its appearance and consistently estimate its value.This is of paramount importance when it comes to comparing numerical simulations to observations.In the former, the Reynolds numbers typical for the latter are not attainable.The present method allows to extract the value of the quadratic drag coefficient from simulations at Reynolds numbers where it is largely dominated by linear friction.An extrapolation to the drag at much higher Reynolds number is then possible.
The drag coefficient obtained ≈1.5×10 −4 (see Table 2) is on the lower end of classical values over smooth surfaces (see Stull, 1988;Schlichting and Gertsen, 2000).The spectral method used (Wirth, 2004) imposes a smoothness of the surface up to the machine precision (10 −16 m).This comes at no surprise as only part of the turbulent small scale motion, responsible for the turbulent viscosity, is explicitly resolved in the two-dimensional non-hydrostatic model, due to the strong anisotropy in the numerical grid (see Sect. 2.3).Drag coefficients over a rough ocean floor are commonly an order of magnitude higher.Current 3-D simulations of only a small part of the a gravity current, but with an isotropic grid resolving the essential turbulent motion, give an increased friction coefficient (preliminary results).
The experiments were performed for moderate surface Rossby numbers as: (i) the numerical resolution in the Navier-Stokes model does not allow for an explicit representation of the fully turbulent Ekman layer (see Sect. 2.3) and (ii) with increasing Rossby number non-hydrostatic effects become important, neither explicitly resolved nor implicitly included (parametrised) in the shallow water model used for the assimilation.
Comparison of free runs of the shallow-water model with the estimated parameter values clearly show that the shallowwater model performs well in the vein, whereas it deteriorates in the friction layer with increasing temperature anomaly.In these areas the local Froud number is larger than unity and the oscillations of the interface and its breaking starts to develop.These non-hydrostatic processes are neither explicitly included in the (hydrostatic) shallow-water model nor are they parametrised.The parametrisation of the dynamics at the interface, including the mixing, entrainment and detrainment, is subject to current research, using higher resolution three dimensional non-hydrostatic simulations of only a small part of the gravity current.Furthermore, to represent the bottom friction in the shallow water model I supposed that unaccelerated linear Ekman layer theory applies to leading order.This assumption is questionable when nonlinear effects, leading to a drag law, become dominant.
I demonstrated that the linear dynamics and the nonlinear departure from it are well captured in the present work.I also showed that data assimilation is a powerful tool in systematically connection models in a model hierarchy.
The parameters of the two most observed and employed friction laws, Rayleigh friction and the drag law, are estimated here.The methodology is not restricted to this two laws, this linear and quadratic laws can also be seen, pragmatically, as the beginning of a Taylor series.

Fig. 11 .
Fig. 11.Cross section of a gravity current with an average geostrophic velocity into the plane.The Coriolis force Fc and the buoyancy force Fg are shown.The figure establishes the force balance between the projection of the Coriolis force F ′c and the projection of the buoyancy force F ′ g onto the topographic slop, for a gravity current on a inclined plane, of angle α, when dissipative processes are neglected.

Fig. 1 .
Fig. 1.Cross section of a gravity current with an average geostrophic velocity into the plane.The Coriolis force F c and the buoyancy force F g are shown.The figure establishes the force balance between the projection of the Coriolis force F c and the projection of the buoyancy force F g onto the topographic slop, for a gravity current on a inclined plane, of angle α, when dissipative processes are neglected.

Fig. 12 .Fig. 2 .
Fig. 12.The drag coefficient cD (straight line) and the effective drag coefficient cD (dashed line), are presented as a function of the Reynolds number (based on the laminar Ekman layer thickness) ReEk.Fig.2. The drag coefficient c D (straight line) and the effective drag coefficient cD (dashed line), are presented as a function of the Reynolds number (based on the laminar Ekman layer thickness) Re Ek .

Table 1 .
Estimated parameter values in the numerical experiments.