Jun-Ichi Yano

Marta Waclwczyk
email: marta.waclawczyk@fuw.edu.pl

Sergej Zilitinkevich
email: sergej.zilitinkevich@fmi.fi

Larry Mahrt
email: mahrt@nwra.com

Danijel Danijel

J.-I Yano

M Wac

Boundary-Layer Meteorology

Keywords: Nondimensionalization scale, Obukhov length, Stably-stratified turbulence, Vertical scale Nondimensionalization scale, Obukhov length, Stably-stratified turbulence, Vertical scale

Attempt is made to derive the Obukhov length as a nondimensionalization scale of partial differential equations governing a turbulent system. When the Richardson number, Ri, is the order of unity or less, this length can be derived as a vertical scale for the nondimensionalizeation based on a balance between the vertical buoyancy flux and the shear-generation terms in the turbulent kinetic energy (TKE) equation. On the other hand, when the Richardson number, Ri, is much larger than unity, the vertical scale is controlled by the buoyancy rather than the background shear. In this case, a characteristic scale is defined by a ratio between the vertical buoyancy eddy flux and a mean vertical advection rate of buoyancy controlled by background stratification. This scale is akin to the buoyancy scale or the external static-stability scale considered in the literature. Moreover, a third characteristic scale is identified as a ratio between the TKE vertical flux and the vertical buoyancy flux. This scale is valid independent of the Richardson number, thus it may be considered an alternative length scale for a turbulence similarity theory.

Introduction

Theories of atmospheric boundary-layer turbulence have been developed by heavily relying on the so-called dimensional analyses [START_REF] Barenblatt | Scaling, Self-Similarity, and Intermediate Asymptotics[END_REF]. This methodology is alternatively called the scaling approach in atmospheric boundarylayer studies, as reviewed by e.g., [START_REF] Holtslag | Scaling the atmospheric boundary layer[END_REF]Nieuwstadt (1986), Foken (2006). Some key variables controlling a given regime of turbulence are first identified, then various characteristic scales of the system (length, velocity, temperature, etc) are determined from these key controlling variables by a dimensional consistency. For example, the Obukhov length [START_REF] Obukhov | Turbulence in an atmosphere with a non-uniform temperature[END_REF]) is defined from the frictional velocity and the buoyancy flux from a dimensional analysis.

For performing a dimensional analysis correctly, a certain ingenuity is required for choosing proper controlling variables of a given system. No systematic methodology exists for choosing them, but the choice is solely based on physical intuitions. A wrong choice of controlling variables can lead to totally meaningless results (cf., [START_REF] Batchelor | Heat convection and buoyancy effects in fluids[END_REF]. With absence of an analytical solution to turbulent flows as well as difficulties in observations and numerical modelling, usefulness of those proposed scales is often hard to judge. As a result, the Obukhov length is hardly a unique choice. There are various efforts to introduce alternative scales as reviewed throughout the text, but the most notably one is a theory based on gradient-based scales [START_REF] Sorbjan | Local structure of stably stratified boundary layer[END_REF][START_REF] Sorbjan | Gradient-based scales and similarity laws in the stable boundary layer[END_REF][START_REF] Sorbjan | Similarity scaling system for stably stratified turbulent flows[END_REF].

However, the dimensional analysis is not a sole possibility of defining characteristic scales of a system. In atmospheric large-scale dynamics, they are typically derived by a nondimensionalization. Arguably, this procedure is more straightforward and formal: we simply introduce characteristic scales for all the variables for nondimensionalzing them. These scales cannot be arbitrary, because we expect that terms in the equation balance each other, thus their orders of magnitudes must match each other. These conditions, in turn, constrain these characteristic scales in a natural manner. Advantage of the nondimensionalization approach is that these scales are defined not only by dimensional consistencies, but also by requirements of balance between the terms in a system. The latter are stronger constraints than the former. The Rossby radius of deformation is a classical example of a characteristic scale identified by a nondimensionalization. This scale characterizes the quasi-geostrophic system (Sec. 3.12, Pedlosky 1987). The depth of the Ekman layer is another such example (Sec. 4.3, Pedlosky 1987).

The present paper attempts to derive the Obukhov length by following this principle of nondimensionalization. For this goal, the principle is applied to a turbulent system. This scale is a core of the celebrated Monin-Obukhov similarity theory [START_REF] Monin | Basic laws of turbulent mixing in the surface layer of the atmosphere[END_REF]. In spite of its importance in describing boundary-layer turbulence, the basis of this Obukhov scaling is often questioned. Especially, various studies suggest that this similarity theory breaks down in a strongly-stratified limit (cf., [START_REF] King | Some measurements of turbulence over an Antarctic ice shelf[END_REF][START_REF] Howell | Surface-layer fluxes in stable conditions[END_REF][START_REF] Mahrt | Stratified atmospheric boundary layers[END_REF]. Various efforts for generalization of the Monin-Obukhov theory already exist (e.g., [START_REF] Zilitinkevich | An extended similarity-theory for the stably stratified atmospheric surface layer[END_REF][START_REF] Zilitinkevich | Third-order transport due to internal waves and nonlocal turbulence in the stably stratified surface layer[END_REF], Zilitinkevich and Easu 2007). However, as far as the authors are aware of, all these efforts are under frameworks of the dimensional analysis. The present paper is going to suggest a procedure beyond those efforts by analysing governing equations of turbulence more directly.

A similar analysis is already performed by [START_REF] Mahrt | Momentum budget of gravity flows[END_REF] more specifically focusing on gravity-wave currents in the boundary layer. Unfortunately, he has only considered a first of step of the nondimensionalization, called "scale analysis" but without actually writing down a nondimensional set of equations.

Advantage of his study is, in turn, in actually solving a closed set of equations.

In contrast, the present study performs the nodimensionalization in a more systematic manner by taking equations for turbulence statistics. However, the chosen set of equations chosen is hardly self-contained. A full nondimensionalization of a turbulence model is performed by Nieuwstadt (1984). However, an adopted model already contains a closure, thus the final nondimensionalization result also depends on the closure. In the present study, we consider turbulence equations without closure for this reason. The present analysis is close in the spirit of [START_REF] George | A similarity theory for the turbulent plane wall jet without external stream[END_REF] in seeking to identify characteristic scales of a given system by directly examining a balance condition in a governing equation, although the latter study does not go through a path of nondimensionalization. Link of the present study to theirs further suggests possibilities of applying various methodologies of multiscale asymptotic expansions to the atmospheric boundary-layer problems, though such attempts are left for future studies.

The paper proceeds as follows. In the next section, the Richardson number is introduced by a nondimensionalization of a linear perturbation problem under a presence of a background vertical wind shear. This example serves for double purposes. First, since the Richardson number plays a key role in the Monin-Obukhov theory, its significance must be better identified in a selfcontained manner (cf., [START_REF] Lobocki | Analysis of vertical turbulent heat flux limit in stable conditions with a local equilibrium, turbulence closure model[END_REF]. Second, it provides a simple demonstration of how characteristic scales can be identified by a nondimensionalization procedure. Here, the analysis can also be performed for a fully nonlinear case as well, but without any major changes in conclusions, except for a consequence of absolute magnitudes of dependent variables can be estimated as a result. This is considered a standard problem in instabilities theories of fluid mechanics, originating from [START_REF] Miles | On the stability of heterogeneous shear flows[END_REF] and [START_REF] Howard | Notes on a paper by John W. Miles[END_REF]. In standard analyses, the Richardson number is introduced in retrospect only after solving the instability problem in dimensional form, with an exception of Sec. 8.1 of [START_REF] Townsend | The Structure of Turbulent Shear Flow[END_REF], which outlines a nondimensinalization procedure for this problem. A nondimensionalization is performed originally in this section, but merely intended to serve as a pedagogical purpose for a demonstration.

The main analyses are presented in Secs. 3 and4, where the turbulentkinetic energy (TKE) and the eddy buoyancy equations are examined, respectively. In these two sections, the nondimensionalization procedure is tailored towards a style of the similarity theory for the turbulence. A full nondimensionalization of this system under a standard procedure is provided in Appendix separately as a reference. Link to the similarity theory is discussed in Sec. 5, and the paper is concluded by final remarks in Sec. 6.

Throughout the paper, Boussinesq approximation is adopted for the analysis for a simplicity. Though an extension of the analysis under an anelastic approximation is relatively straightforward, there is no much advantage in this generalization, but only with a consequence of making the analysis more involved. For the same reason, a two-dimensional system is considered in Sec. 2 with x and z taken as horizontal and vertical coordinates, respectively.

The bar, ¯, and the prime, ′, signs are used throughout for designating the mean and the deviation from the former. The former is assumed to depend only on height, z. The latter corresponds to the turbulence fluctuation when a fully nonlinear problem is considered, as in Secs. 3 and 4. On the other hand, under a linear stability analysis in Sec. 2, the primed quantities designate the perturbation variables. Due to the linearization, the perturbation variables may grow to infinity under an unstable situation, whereas the turbulent fluctuations are bounded by full nonlinearity. Nevertheless, the formal definition of the prime sign itself does not change over these sections.

A Linear Perturbation Problem and the Richardson Number

As already stated, only a perturbation problem is considered in this section by neglecting the nonlinearity by following a standard shear instability problem [START_REF] Miles | On the stability of heterogeneous shear flows[END_REF][START_REF] Howard | Notes on a paper by John W. Miles[END_REF]. Also by following it, the viscosity is neglected in this section. Thus, the governing equations of the problem are given by:

∂u ′ ∂t + w ′ dū dz + ū ∂ ∂x u ′ = - ∂φ ′ ∂x , (1a)
∂w ′ ∂t + ū ∂ ∂x w ′ = - ∂φ ′ ∂z + b ′ , ( 1b
)
∂b ′ ∂t + w ′ d b dz = 0, ( 1c
)
∂u ′ ∂x + ∂w ′ ∂z = 0. ( 1d
)
Here, u and w are the horizontal and the vertical components of the velocity, φ the pressure divided by a reference density, b the buoyancy, which may be evaluated from the potential temperature, θ, as b = gθ/θ 0 , where g is the acceleration of the gravity, θ 0 is a reference value of the potential temperature.

For the moist atmosphere, the potential temperature must be furthermore replaced by the virtual potential temperature for accounting for a contribution of moisture to the buoyancy.

Nondimensionalizations are performed on the variables by designating the nondimensional variables by the dagger, †. The nondimensionalization scales for the dependent variables are designated by the subscript * . Thus,

u ′ = u * u † , w ′ = u * w † , φ ′ = φ * φ † , b ′ = b * b † , ∂ ∂t = 1 t * ∂ ∂t † , ∂ ∂x = 1 z * ∂ ∂x † , ∂ ∂z = 1 z * ∂ ∂z † ,
and also

ū = ū * ū † , dū dz = dū dz * dū † dz † , d b dz = d b dz * d b † dz † .
Here, the velocity and the spatial scales are not differentiated in horizontal and vertical directions for maintaining the order of magnitude of the pressure gradient the same in both directions, thus x * = z * and w * = u * . Also note that the background-state gradients, dū/dz and d b/dz, are nondimensionalized by directly using the scales for these gradients, rather than by ū/z * and b/z * , respectively, for the purpose deriving a standard definition of the Richardson number below.

By substitutions of those expressions into the equations, we find the following relations for the nondimensionalization scales:

t * = z * ū * , φ * = ū * u * , b * = ū * u * z * . (2a, b, c)
Here, the last relation can be interpreted as a definition of the length scale by re-writing it as:

z * = ū * u * b * , (3)
if all the variables in the right hand side can be considered to be prescribed.

As a result, a nondimensionalized set of equations is given by

∂u † ∂t † + z * ū * dū dz * dū † dz † w † + ū † ∂u † ∂x † = - ∂φ † ∂x † (4a) ∂w † ∂t † + ū † ∂w † ∂x † = - ∂φ † ∂z † + b † (4b) ∂b † ∂t † + z * ū * 2 d b dz * d b † dz † w † = 0 (4c) ∂u † ∂x † + ∂w † ∂z † = 0 (4d)
Note that by choosing the nondimensionalization scales by Eqs. (2a, b, c), there is no constant factor in front of almost all the terms in the equations, except for the two in Eqs. (4a) and (4c). Here, by choosing a length scale, z * , in an appropriate manner, we can remove a constant factor from one of these two terms, but not from both of them.

There are two options for choosing the length scale, z * : (i) by setting a factor in front of the second term of Eq. (4a) unity, or (ii) by setting a factor in front of the second term in Eq. ( 4c) unity. The option (i) amounts to set the spatial scale (shear scale) as that of the background wind shear, and we obtain

z * = z * u ≡ ū * /(dū/dz) * . (5a)
The option (ii) leads to

z * = z * b ≡ ū * /(d b/dz) 1/2 * . (5b)
We may call the latter the buoyancy-gradient scale.

As a result, the set of equations also reduces to with the option (i):

∂u † ∂t † + d ū † dz † w † + ū † ∂u † ∂x † = - ∂φ † ∂x † , ∂w † ∂t † + ū † ∂w † ∂x † = - ∂φ † ∂z † + b † , ∂b † ∂t † + Ri d b † dz † w † = 0, ∂u † ∂x † + ∂w † ∂z † = 0,
and with the option (ii):

∂u † ∂t † + Ri -1/2 dū † dz † w † + ū † ∂u † ∂x † = - ∂φ † ∂x † , ∂w † ∂t † + ū † ∂w † ∂x † = - ∂φ † ∂z † + b † , ∂b † ∂t † + d b † dz † w † = 0, ∂u † ∂x † + ∂w † ∂z † = 0.
Here, Ri is the Richardson number defined by a ratio of the two characteristic scales:

Ri = z * u z * b 2 = (d b/dz) * (dū/dz) 2 * . ( 6
)
We find that when the shear is more dominant than the buoyancy gradient (stratification), i.e., Ri < 1, the scaling based on the shear scale, z * u , (Eq. 5a) is relevant, and when the stratification is more dominant than the shear, i.e., Ri > 1, the scaling based on the buoyancy-gradient scale, z * b , (Eq. 5b) becomes relevant. We expect that these Ri-dependent characteristics of the system are still valid also for fully turbulent regimes, that are going to be addressed in the next two sections. An equivalent nondimensionalization, as in this section, is performed, separately, in the Appendix for a full turbulence system considered in Secs. 3 and 4, and essentially the same conclusion is drawn.

In this manner, we have demonstrated how naturally characteristic scales (not only the length scale) of a system can be determined by a nondimensionalization. A question that we are going to address in the next two sections is whether the Obukhov length can be derived in a similar manner for fully turbulent flows.

Furthermore, the result obtained in this section already has implications in the boundary-layer turbulence, because the identified characteristic scales, z u * and z b * , (Eqs. 5a, b) are expected to characterize the typical size of eddies of given regimes, and thus, also characterize the resulting mixing lengths, l * .

In this respect, it may be worthwhile to note that, for example, [START_REF] Grisogono | Generalizing 'z-less' mixing length for stable boundary layers[END_REF] propose to use two different mixing lengths,

l * = v ′ * (dū/dz) * , (7a)
l * = v ′ * (d b/dz) 1/2 * , (7b)
depending on the Richardson number, Ri. Here, v ′ * is a scale for the velocity.

A similar scale (buoyancy scale), defined by

l * = w ′ * (d b/dz) 1/2 * , (8)
is introduced by [START_REF] Stull | Inversion rise model based on penetrative convection[END_REF], [START_REF] Zeman | Parameterization of the turbulent energy budget at the top of the daytime atmospheric boundary layer[END_REF], [START_REF] Brost | A model study of the stably stratified planetary boundary layer[END_REF]. [START_REF] Hunt | Some observational structure in stable layers[END_REF], in turn, suggest by field data analysis that the buoyancy scale (Eq. 8) characterizes both the vertical heat transport and the temperature-variance production in the stably-stratified boundary layer.

These definitions reduces to z * u and z * b , respectively, with small and large

Richardson numbers by re-setting as v ′ * = ū * and w ′ * = ū * . This condition is expected to be satisfied when a system is fully turbulent.

3 Turbulent System: Turbulent Kinetic Energy Equation

Obukhov Length

In the following two sections, nondimensionalization of stably-stratified turbulence system is considered. The goal is to derive the Obukhov length as a natural consequence of nondimensionalization, in a similar manner as demonstrated in the last section how a characteristic length of a system is identified.

Thus, we focus on the equations containing the vertical buoyancy flux, w ′ b ′ , and the vertical momentum stress, u ′ w ′ . The Obukhov length is defined by a ratio of fractional powers of those two quantities:

L = (u ′ w ′ ) 3/2 * (w ′ b ′ ) * . (9)
Here, (u ′ w ′ ) * and (w ′ b ′ ) * are the characteristic scales for the vertical momentum stress and the vertical buoyancy flux, respectively. Note that unlike in the dimensional analysis, under a nondimensionalization procedure, these scales do not necessarily refer to actual values at a particular vertical level (say, the surface). Even when such a choice to be made, the chosen values must also be representative for the whole vertical stretch of the boundary layer.

This definition is often further simplified into

L = u * 3 (w ′ b ′ ) *
by introducing a frictional velocity, u * , defined by

u * 2 = (u ′ w ′ ) * . ( 10
)
Here, note that in a standard nondimensionalization procedure, the orders of magnitudes of those flux terms are estimated by

(u ′ w ′ ) * = u * w * , (w ′ b ′ ) * = w * b * .
By substituting these two expressions into (9),

L = u 3/2 * w 1/2 * b * .
If an isotropic scaling (i.e., w * = u * ) can further be assumed as in the last section, and also ū * = u * , the last expression reduces to Eq. ( 3) with L = z * .

In this manner, taking the Obukhov length as an example here, we see that the turbulence length scales assume more than what are typically assumed for the nondimensionalization scales. See Secs. 3.7 and 5 for the further discussions.

Turbulent Kinetic Energy Equation

In this section, we consider the turbulent kinetic energy (TKE) equation, because it contains both the vertical buoyancy flux, w ′ b ′ , and the vertical momentum stress, u ′ w ′ . Based on a result from the last section, we expect that a closed scaling for the nondimensionalization is obtained by considering this equation when the shear effect is strong (i.e., Ri ≤ 1). On the other hand, when the stratification becomes strong enough (i.e., Ri ≫ 1), we expect that we also have to consider the eddy-buoyancy equation, which also contains a vertical buoyancy flux term. The latter will be considered in the next section.

By following a standard formulation of the atmospheric boundary-layer turbulence, only the vertical flux terms are considered assuming horizontal homogeneity. This is solely for simplifying the analysis focusing on the goal of deriving the Obukhov length. Horizontal heterogeneity is expected to be important for some stable-stratified atmospheric turbulent flows, but this extension is left for a future study.

A standard TKE equation (e.g., Deardorff 1983) is given by:

∂ ∂t v ′2 2 = w ′ b ′ -u ′ w ′ ∂ ū ∂z - ∂ ∂z w ′ (v ′2 + φ ′ ) -ε (11)
Here, the overbar designates a horizontal average, v is a velocity vector, and ε the dissipation rate.

3.3 The Balance:

w ′ b ′ ∼ (u ′ w ′ )dū/dz
We first consider a balance between two terms that involve the vertical buoyancy flux, w ′ b ′ , and the vertical momentum stress, u ′ w ′ , respectively, because this balance is likely to lead to a derivation of the Obukhov lenght. This amounts to consider a balance between the first and the second terms. This balance requirement leads to:

(w ′ b ′ ) * = (u ′ w ′ ) * dū dz * . ( 12
)
This condition may be used to estimate a vertical scale, z * , by introducing a typical change, ū * , of the background wind over this scale so that

dū dz * = ū * z * . ( 13
)
By substituting this relation back to the balance condition ( 12), we obtain an estimate of the length scale given by

z * = ū * (u ′ w ′ ) * (w ′ b ′ ) * . ( 14
)
Our next goal is to try to show the equivalence of this length scale with that of Obukhov.

The

Balance: (u ′ w ′ )dū/dz ∼ ∂(w ′ v ′2 )/∂z
This length scale (14) reduces to the Obukhov length (9) if we can set ū * = (u ′ w ′ )

1/2 * .

To get an answer to this question, we need an estimate of the windshear strength, ū * , relative to the eddy. Thus, this amounts to consider a balance between the second and the third terms.

This balance condition is given by

(u ′ w ′ ) * dū dz * = (w ′ v ′2 ) * z * .
To obtain a more explicit expression for (dū/dz) * in terms of the turbulence fluctuations, we further note

(u ′ w ′ ) * = u * w * , (15a)
(w ′ v ′2 ) * = u 2 * w * . (15b)
By substituting Eqs. (15a, b) into the balance condition above, it reduces to

dū dz * = u * z * , (16)
and ū * = u * . Thus,

z * = u * (u ′ w ′ ) * (w ′ b ′ ) * . ( 17
)
The stress term can further be re-written as

(u ′ w ′ ) = ǫu 2 *
by introducing an aspect ratio, ǫ, of the system. As a result,

u * = ǫ-1/2 (u ′ w ′ ) 1/2 . ( 18
)
By further substituting ( 18) into (17), we finally obtain

z * = ǫ-1/2 (u ′ w ′ ) 3/2 * (w ′ b ′ ) * . ( 19
)
If the turbulence is isotropic, the length scale ( 19) reduces to the Obukhov length (9). However, keep in mind that stably-stratified turbulence is often observed to be quasi two dimensional (i.e., ǫ ≪ 1), thus the Obukhov length may underestimate the vertical scale.

Most importantly, keep in mind that the balance ( 16) is possible only if the wind-shear is strong enough, and the Richardson number, Ri, is the order of unity or less. When the Richardson number, Ri, is very small, this balance is no longer valid (the wind-shear term drops as O(Ri)), and the perturbationbuoyancy equation must be considered instead, as in Sec. 4, in defining the length scale of the system.

3.5 The Balance:

w ′ b ′ ∼ ∂(w ′ v ′2 )/∂z
In completing the analysis, we also consider the balance between the first and the third terms, because this term provides an alternative estimate of the vertical scale, z * . It is given by

z * = (w ′ u ′2 ) * (w ′ b ′ ) * . ( 20
)
Note that the estimate of the length scale by Eq. ( 20) is valid regardless of the magnitude of the Richardson number, Ri, of the system, unlike Eq. ( 19).

Deardorff Velocity Scale

Eq. ( 20) can further be re-written by noting that

(w ′ u ′2 ) * = u 2 * w * = ǫu 3 * = ǫ-2 w 3 *
Rearrangement after the substitution leads to a velocity scale

w * = [ǫ 2 z * (w ′ b ′ ) * ] 1/3 (21)
This velocity scale reduces to one introduced by [START_REF] Deardorff | Convective velocity and temperature scales for the unstable planetary boundary layer and for Rayleigh convection[END_REF] and [START_REF] Tennekes | Free convection in the turbulent Ekman layer of the atmosphere[END_REF], when we set ǫ = 1 and z * the boundary-layer depth. However, it may be useful to realize that the same expressions can also be derived by adopting eddy-diffusion formulations. Under the latter, the left hand sides in Eqs. (15a, b) may be represented as

(u ′ w ′ ) * = -ν t dū dz (22a) (w ′ v ′2 ) * = -ν t dū 2 dz (22b)
Here, the eddy-diffusion coefficient, ν t , may further be written as

ν t = u * l *
in terms of a velocity scale, u * , and a mixing length, l * . By substituting this expression into Eqs. (22a, b), we finally obtain

(u ′ w ′ ) * = u * l * ū * z * f † (z † ) (23a) (w ′ v ′2 ) * = u * l * ū2 * z * g † (z † ) (23b)
Here, f † (z † ) and g † (z † ) are universal functions (cf., Sec. 5 for further discussions). By normalizing them as, say, f † (1) = 1 and g † (1) = 1, and also However, here, more steps are required to reach the same conclusion.

Turbulent System: Buoyancy Perturbation Equation

The analysis in Sec. 2 suggests that as the Richardson number, Ri, becomes larger than unity, the shear term in the TKE equation (or the momentum equation) becomes less important, so that it loses a basis of using this term for deriving the Obukhov length from Eq. ( 14). As the Richardson number increases, in turn, the buoyancy stratification, d b/dz, plays a more important role. For considering the contribution of this term, in this section, we consider the buoyancy perturbation equation:

∂b ′ ∂t + w ′ d b dz + ∂w ′ b ′ ∂z = Q ′ . (24)
Here, Q ′ is a perturbation diabatic heating.

By considering a balance between the second and the third terms in the left hand side, we obtain as a length scale:

z * = (w ′ b ′ ) * w * (d b/dz) * . ( 25
)
This is expected to be a characteristic scale of a turbulent flow when the perturbation buoyancy becomes a dominant under a strong stratification, being consistent with a more general analysis in the Appendix.

This scale is somehow akin to the external static-stability scale

L N = (u ′ w ′ ) 1/2 * (d b/dz) 1/2 * .
as introduced by [START_REF] Kitaigorodskii | A note on similarity theory for atmospheric boundary layers in the presence of background stable stratification[END_REF], and considered, especially, by [START_REF] Zilitinkevich | Resistance and heat-transfer for stable and neutral planetary boundary layers: Old theory advanced and reevaluated[END_REF]. More specifically, [START_REF] Zilitinkevich | An extended similarity-theory for the stably stratified atmospheric surface layer[END_REF] suggests a nondimensional parameter, L/L N , to define a transition from a regime dominated by the Obukhov lenght, L, to L N . A link to the buoyancy-gradient scale introduced by Eq. ( 5b) is also noted. Recall that the latter is further linked to the buoyancy scale introduced by Eqs. (7b, c). [START_REF] Sorbjan | Local structure of stably stratified boundary layer[END_REF][START_REF] Sorbjan | Gradient-based scales and similarity laws in the stable boundary layer[END_REF][START_REF] Sorbjan | Similarity scaling system for stably stratified turbulent flows[END_REF], in turn, develops his gradient-based similarity theory based on the buoyancy scale. [START_REF] Zilitinkevich | Resistance and heat-transfer for stable and neutral planetary boundary layers: Old theory advanced and reevaluated[END_REF] argue that the scale, L N , becomes relevant when the vertical eddy heat flux is small. In turn, we argue that the scale defined by Eq. ( 25) becomes relevant when the Richardson number, Ri, is large enough. The scale, L N , can be obtained by assuming that Ozmidov scale,

L o = ε 1/2 (d b/dz) 3/2 * , is equal to another length scale, z * = (u ′ w ′ ) 3/2 * /ε,
which defines the characteristic turbulence length [START_REF] Grisogono | Generalizing 'z-less' mixing length for stable boundary layers[END_REF].

It transpires that the scale, L N , is purely based on a dimensional analysis.

In contrast, the scale ( 25) is derived from an actual balance in an equation.

These two scales become equivalent when all the flux scales are re-written in terms of the scales of more basic variables as in Eqs. (15a, b), and also the relation (3) is invoked. We expect that values of these two scales are also similar numerically in practice. Note furhter that under this equivalence, the nondimensional parameter, L/L N , also reduces to the Richardson number, Ri.

However, as emphasized in the next section, in employing a length scale in a context of a dimensional analysis leading to a similarity theory, it must be exact numerical value rather than just an order-of-magnitude estimate, as in a standard nondimensionalization.

Discussions: Link to the Similarity Theory

A basis for the similarity theory may be provided from the nondimensionalization of the system considered in the last two sections in the following manner. As an example, let us consider the TKE equation ( 11). Note that the TKE equation is hardly self-contained. However, for the sake of a heuristic argument, let us suppose that every variable in Eq. ( 11) can be determined self-consistently by solving it. Alternatively, we may just suppose, though not feasible in reality, that all the necessary higher-order moments are also nondimensionalized in a similar manner, and all the necessary equations are solved to obtain all the variables in Eq. ( 11).

Steady solutions for these nondimensionalized variables are given in terms of nondimensional functions, say, f † (z † ), g † (z † ), q † (z † ), etc. as

w †′ b †′ = f † (z † ) (26a) u †′ w †′ = g † (z † ) (26b) ∂ ū † ∂z † = q † (z † ), (26c)
etc. After dimensionalizations,

w ′ b ′ = (w ′ b ′ ) * f † ( z z * ) (27a) u ′ w ′ = (u ′ w ′ ) * g † ( z z * ) (27b) ∂ ū ∂z = ∂ ū ∂z * q † ( z z * ) (27c)
Furthermore, let us suppose that the only necessary boundary conditions required for determining these variables are their surface values. In that case, As a result, general solutions to the system are given by

we set (w ′ b ′ ) * = (w ′ b ′ ) 0 , (u ′ w ′ ) * = (u ′ w ′ ) 0 ,
w ′ b ′ = (w ′ b ′ ) 0 f † ( z z * ) (28a) u ′ w ′ = (u ′ w ′ ) 0 g † ( z z * ) (28b) ∂ ū ∂z = ∂ ū ∂z 0 q † ( z z * ) (28c)
Based on the arguments so far, we may conclude that the nondimensional functions, f † (z/z * ), g † (z/z * ), q † (z/z * ), are universal only depending on the nondimensionalization scale, z * . This is an essence of the similarity theory (cf., [START_REF] Sorbjan | Structure of the Atmospheric Boundary Layer[END_REF].

Here, for this last statement to be valid in strict manner, it is not sufficient that the nondimensionalization scale, z * , simply measures a characteristic scale of a system, merely as an estimate of an order of magnitude, as the nondimensionalization intends to do. Instead, this length scale, z * , must be re-scaled in a precise manner as the system environment changes so that the universality of the functions f , g † , and q † is maintained.

Being consistent with this requirement for developing a similarity theory, in the last two sections, attempt is made to define the length scale, z * , of the system in a more strict manner based on the flux scales, (u ′ w ′ ) * , (w ′ b ′ ) * etc.

rather than simply in terms of the variable scales, u * , w * , etc. For applying a similarity theory in a strict manner, these flux values must also be defined at a specific vertical level, probably, most conveniently at the surface.

In this respect, under a systematic nondimensionalization procedure, the present study has failed to derive the Obukhov length, L, (Eq. 9) in any strict manner as a nondimensionalization scale, z * . The closest we have obtained is Eq. ( 14), which defines the length scale in terms of both the wind stress magnitude, (u ′ w ′ ) * , and the wind-shear scale, ū * . It reduces to the Obukhov length only under an extra scale argument (15a, b) for isotropic turbulence (i.e., ǫ = 1). We should realize that this reduction is hardly exact from a point of view of the similarity theory, either, thus in applying this scale to the universal functions, f † , g † , and q † , some minor adjustment factors must be multiplied on the nondimensionalized vertical coordinate, z † .

Most importantly, this scaling is valid only when the Richardson number is Ri ≤ 1. When Ri ≫ 1, with strong stratification, the buoyancy-gradient scale defined by Eq. ( 25) becomes more relevant rather than the shear-based scale ( 14).

However, the TKE equation identifies an alternative vertical scale defined by Eq. ( 20), which is expected to be valid independent of the value of the Richardson number, thus also independent of the turbulence regimes. This new length scale is not remotely different from the Obukhov length, but obtained by replacing the frictional velocity, u * , by the scale, (w ′ u ′2 )

1/3 * , for the vertical flux of the turbulent kinetic energy.

Summary and Further Remarks

The present study has derived the Obukhov length as a result of a balance between the vertical buoyancy flux and the shear-generation terms in the TKE equation. However, this balance is valid only with the Richardson number, Ri, of the order of unity or less. When the Richardson number is much larger than unity, the system is characterized by another length scale defined by a ratio between the vertical buoyancy eddy flux and a vertical advection rate of buoyancy controlled by background stratification (Eq. 25). This scale is akin to the buoyancy scale introduced by [START_REF] Stull | Inversion rise model based on penetrative convection[END_REF], [START_REF] Zeman | Parameterization of the turbulent energy budget at the top of the daytime atmospheric boundary layer[END_REF], [START_REF] Brost | A model study of the stably stratified planetary boundary layer[END_REF], and the external static-stability scale considered by [START_REF] Zilitinkevich | Resistance and heat-transfer for stable and neutral planetary boundary layers: Old theory advanced and reevaluated[END_REF].

The present nondimensionalization analysis has also identified a new length scale defined as a ratio between the the TKE vertical flux, (w ′ u ′2 ) * , and the vertical buoyancy flux, (w ′ b ′ ) * (Eq. 20). This scale may be adopted as an alternative length scale for developing a boundary-layer turbulence similarity theory. A next step would be to examine this possibility based on field-campaign data sets.

An important aim behind the present study has been to demonstrate how a characteristic scale of a turbulent system can be identified directly by nondimensionalization of a partial differential-equation system describing turbulence. For this demonstrative purpose, the analysis has been performed with the simplest possible turbulence system, assuming a horizontal homogeneity.

Also only a limited set of equations is examined, namely, the TKE and the buoyancy-perturbation equations. The choice is made specifically with a goal in mind of identifying the Obukhov length so that equations contain the dimensional parameters used in its definition.

An equivalent nondimensionalization analysis with a full turbulence equation system is still to be performed. Nevertheless, the present preliminary analysis already suggests a fruitfulness of such an investigation. It is expected that different turbulent regimes are identified by changing orders of magnitudes of the Richardson number, as the present study has already suggested.

Such an analysis is expected to provide a more solid theoretical basis for interpreting the various different turbulent regimes phenomenologically identified for the stably-stratified turbulence (cf., [START_REF] Mahrt | Stratified atmospheric boundary layers[END_REF]. Various further generalizations are equally feasible. The present analysis has been performed under an assumption of quasi-stationarity of the system. However, some of the turbulence regimes under stable stratification may be fundamentally transient (cf., [START_REF] Caughey | Turbulence in the evolving stable boundary layer[END_REF]. A role of horizontal heterogeneity is another aspect to be investigated, especially in a context of stably-stratified turbulence. For example, under certain situations, the horizontal heat transport, a term that is often neglected in theoretical studies, becomes a key process in heat budget (e.g., [START_REF] Wittch | The nocturnal boundary layer over the Northern Germany: An observational study[END_REF]. A role of anisotropy of the flow with ǫ ≪ 1 is still to be carefully examined as well.

Another aspect, that is not discussed herein, is a role of boundary conditions in solving the turbulence problems. When our focus is on a layer close enough to the surface (e.g., surface layer), a contribution from a top of the planetary layer may be neglected, as a basic premise of the Monin-Obukhov theory as well as in subsequent generalizations. However, when a problem concerns a whole depth of the boundary layer, the depth of the boundary layer becomes another parameter to be considered. As pointed out by e.g., [START_REF] Holtslag | Scaling the atmospheric boundary layer[END_REF], the problem must be solved by explicitly taking into account of a condition at the top of the boundary layer.

A critical difference must also be recognized between a typical dimensional analysis performed in the turbulence studies and the nondimensionalization procedure considered in the present study. It is common in turbulence studies to take characteristic scales (not only the length scale, but more generally) to be functions of height, as manifested as a local similarity theory (Nieuwstadt 1984, Sec. 4.4.1, Sorbjan 1989;[START_REF] Mahrt | Stratified atmospheric boundary layers[END_REF]. For example, a frictional velocity can be evaluated as a function of height, under a dimensional analysis, and the given result is used for furthermore evaluating the Obukhov length.

On the other hand, a single fixed scale must be adopted for each variable in nondimensionalization. If the velocity is ever nondimensionalized by a scale that depends on height, a height-dependence of the dimensional velocity must somehow be partitioned between a nondimensional velocity and a velocity scale for a nondimensionalization. Such a partitioning is hardly trivial, thus usually not considered.

A major exception to this rule is when variables are defined purely locally, without influence of their spatial gradients (e.g., with absence of transport processes). In this case, nondimensionalization may be performed independently at every spatial point. Conversely, it suggests that a local similarity theory is valid in a strict sense, only if a system is locally in equilibrium at each vertical level without significant effects of vertical transport, as theoretically expected as a z-less stratification state (Neustadt 1984). In other words, a local similarity theory is strictly valid only for the flows with strong two dimensionality.

The last remark further suggests how a validity of a given similarity theory can be verified by a more direct nondimensionalization of a given system, more generally speaking. Further investigations along this line are much anticipated.

Acknowledgements Comments and encouragements by Szymon Malonowski are must appreciated.

Appendix: Full Nondimensionalization of the Turbulent System

The appendix presents a full nondimensionalization of the system consisting of Eqs. ( 11) and (24). By following the same procedure as in Sec. 2, we find the following relationships between the nondimensionalization scales: In analogous manner as in Sec. 2, we obtain a nondimensionalized system:

∂ ∂t † v †′2 2 = w †′ b †′ - z * u * ∂u ∂z * u †′ w †′ ∂ ū † ∂z † - ∂ ∂z † w †′ (v †′2 + φ †′ ) -ε † (30a) ∂b †′ ∂t † + z * u * 2 db dz * w †′ d b † dz † + ∂w † ′ b † ′ ∂z † = Q †′ (30b)
Here, the vertical scale, z * , is left unspecified. As in Sec. 2, we have two

3. 7

 7 Link to Eddy-Diffusion Formulation Those readers who are with a strong background in boundary-layer meteorology, Eqs. (15a, b) may not be immediately clear. Here, the right hand sides are the expressions used in the standard nondimensionalization procedure, as considered in Sec. 2. As already remarked at the end of Sec. 3.1, the dimensional analysis in the boundary-layer meteorology, more specific scale variables, such as fluxes, are considered for deriving the characteristic scales of a system. On the other hand, under a nondimensionalization procedure, only the orders of magnitudes of the variables in concern, thus if u ∼ u * , w ∼ w * , etc, we can immediately write as Eqs. (15a, b). Keep in mind that, by definition, the values in the left-hand sides are also only the orders of magnitude estimates. They do not correspond at all to any actual flux values at any vertical levels, as expected in standard boundary-layer similarity theories.

  assuming l * = z * , ū * = u * , ǫ = 1, Eqs. (23a, b) reduce to Eqs. (15a, b).

  and (∂ ū/∂z) * = (∂ ū/∂z) 0 with the subscript, 0, designating the surface values. Here, those eddy flux values are more precisely defined at the top of the viscous boundary layer, but these can also be equated with the actual values of the fluxes from the surface.

  options to choose: (i) (z * /u * )(∂u/∂z) * = 1, or (ii) (z * /u * ) 2 (db/dz) * = 1. Each leads to an option of (i) the wind-shear scale, z * = u * /(∂u/∂z) * or (ii) the buoyancy-gradient scale, z * = u * /(db/dz) 1/2 * . In respective cases, we obtain (i) (z * /u * ) 2 (db/dz) * = Ri, and (ii) (z * /u * )(∂u/∂z) * = Ri -1/2 . Recall that the Richardson number, Ri, is defined by Eq. (6).


O El Ganaoui
email: omar.elganaoui@irisa.fr

X Morandi

S Duchesne

P Jannin

Preoperative Brain Shift: Study of Three Surgical Cases

Keywords: Neurosurgical Procedures, Brain Shift, MIGNS, Surgical Workflow

In successful brain tumor surgery, the neurosurgeon's objectives are threefold: (1) reach the target, (2) remove it and (3) preserve eloquent tissue surrounding it. Surgical Planning (SP) consists in identifying optimal access route(s) to the target based on anatomical references and constrained by functional areas. Preoperative images are essential input in Multi-modal Image Guided NeuroSurgery systems (MIGNS) and update of these images, with precision and accuracy, is crucial to approach the anatomical reality in the Operating Room (OR). Intraoperative brain deformation has been previously identified by many research groups and related update of preoperative images has also been studied. We present a study of three surgical cases with tumors accompanied with edema and where corticosteroids were administered and monitored during a preoperative stage [t 0 , t 1 = t 0 + 10 days]. In each case we observed a significant change in the Region Of Interest (ROI) and in anatomical references around it. This preoperative brain shift could induce error for localization during intervention (time t S ) if the SP is based on the t 0 preoperative images. We computed volume variation, distance maps based on closest point (CP) for different components of the ROI, and displacement of center of mass (CM) of the ROI. The matching between sets of homologous landmarks from t 0 to t 1 was performed by an expert. The estimation of the landmarks displacement showed significant deformations around the ROI (landmarks shifted with mean of 3.90 ± 0.92 mm and maximum of 5.45 mm for one case resection). The CM of the ROI moved about 6.92 mm for one biopsy. Accordingly, there was a sizable difference between SP based at t 0 vs SP based at t 1 , up to 7.95 mm for localization of reference access in one resection case. When compared to the typical MIGNS system accuracy (2 mm), it is recommended that preoperative images be updated within the interval time [t 1 ,t S [ in order to minimize the error correspondence between the anatomical reality and the preoperative data. This should help maximize the accuracy of registration between the preoperative images and the patient in the OR.

INTRODUCTION

This study is set within the context of brain tumor surgery. In this paper, we demonstrate the necessity to update preoperative images in order to allow the recovery of registration error for Multi-modal Image Guided Neurosurgery (MIGNS) systems due to a phenomenon that we call preoperative brain shift. In successful surgery, the neurosurgeon's objectives are threefold: (1) reach the target, (2) remove it and (3) preserve eloquent tissue surrounding it. Surgical Planning (SP) consists in identifying optimal access route(s) to the target based on anatomical references and constrained by functional areas. [START_REF] Jannin | Integration of sulcal and functional information for multimodal neuronavigation[END_REF][START_REF] Barillot | Image guidance in neurosurgical procedures, the Visages point of view[END_REF] MIGNS systems have been developed and are being used to help reach these objectives. A key component of MIGNS is the accurate registration of preoperative images to the intraoperative coordinate system of the patient. [START_REF] Peters | Image-guidance for surgical procedures[END_REF] The first registration step consists in a rigid transform, based on assumptions that are not verified during neurosurgical procedures. Intraoperative brain deformation has been previously identified by many research groups and related update of preoperative images has also been studied. [START_REF] Miga | Model-updated image guidance: initial clinical experiences with gravity-induced brain deformation[END_REF][START_REF] Roberts | Intraoperative brain shift and deformation: a quantitative analysis of cortical displacement in 28 cases[END_REF] One study [START_REF] Hartkens | Measurement and analysis of brain deformation during neurosurgery[END_REF] has shown that the surface of the brain is deformed by up to 25mm after skull opening, which could lead to substantial errors in MIGNS systems. To our knowledge, there is no study that discuss the possible deformations before skull opening. In particular, our study addresses the problem of possible deformations in the case of tumors associated with edema at the preoperative stage. Thus, in this paper we introduce the concept of preoperative brain shift. Tissue edema is a pathological condition involving an increase of the amount of fluid in tissue. [START_REF] Grasso | Assessment of human brain water content by cerebral bioelectrical impedance analysis: a new technique and its application to cerebral pathological conditions[END_REF] Brain edema represents a disturbance of the volume equilibrium which must be compensated by a reduction in other fluid and blood compartments. When this compensation is inadequate, tissue pressure and intracranial pressure increase. For standard hospital protocol, after the diagnosis, the edema is treated by administration of corticosteroids during the preoperative phase until intervention time and no update of images is requested before the intervention. Let us define t 0 as the time of preoperative image acquisition and t S the time of intervention. During [t 0 ; t S [, the patient is administered corticosteroids to regularize the intracerebral pressure and to improve neurological deficit due to edema. Since no imaging is typically performed during this interval, the MIGNS system registration, with 2 mm 8 of accuracy, is based on preoperative images acquired at time t 0 . Our hypothesis is that there will be deformations of the ROI induced by corticosteroids administration during preoperative stage [t 0 ; t S [. The edema should shrink, while at the same time the tumor will likely grow. This will induce a significant displacement of neighboring structures over a large region of interest; and to illustrate that point, we present a study of three patients with tumors and severe accompanying edema where corticosteroids were administered and monitored during a preoperative stage [t 0 , t 1 = t 0 + 10 days]. Our objective is to capture the shift in closed cranium, in these cases where parameters are controlled (tumor localization, edema volume, tumor volume, and corticosteroids administration) without surgical manipulation and mechanical deformation relied to skull opening. To initialize properly the MIGNS systems by correcting the preoperative brain shift implies proper initialization of the well known brain shift phenomenon.

MATERIAL AND METHODS

Clinical data and images protocol acquisition

The studied cases were three patients with tumors presenting an important surrounding edema: the first patient presented tumor distributed bilaterally with foci in both left and right hemispheres, the second patient revealed right frontal tumor, and the third patient presented ventricular tumor. The treatment consisted in tumor biopsy for the first and third patient, and tumor resection for the second patient. Our imaging protocol included 3T-MR imaging (Philips Medical Systems) using: (1) three-Dimensional (3D) T1-weighted MRI (T R = 9.89 ms, T E = 4.59 ms, 1 voxel = 1 mm 3 ), (2) 3D T1-weighted Gadolinium enhanced (T R = 9.89 ms, T E = 4.59 ms, 1 voxel 1 mm 3 ) to increase the detection of the tumor and (3) 2D axial T2-weighted Flair (T R = 11000 ms, T E = 125 ms, resolution : 0.34 × 0.34 × 4 mm 3 ) to maximize information about the edema. In respect to our hypothesis, we set the time t 1 to 10 days after the first preoperative acquisition (t 1 = t 0 +10 days).

Image Processing WorkFlow (IPWF)

The IPWF was broken down in two parts: (1) segmentation of the structures in the different MRI sequences; and (2) reformating of all structures in the same reference frame (T1-weighted at acquisition time t 0 ) for comparison. The IPWF (Figure 1) included:

1. Image preprocessing: we corrected the images intensity inhomogeneities 9 in order to improve the segmentation process. To minimize the impact of noise on the image processing performance, we denoised the images using NL-means algorithm. [START_REF] Coupé | Fast Non Local Means Denoising for 3D MR Images[END_REF] 2. Image registration: we performed rigid image registration by maximization of the mutual information between the images using the NEWUOA algorithm. [START_REF] Wiest-Daesslé | Evaluation of a new optimisation algorithm for rigid registration of MRI data[END_REF] We note the matrix of rigid-body registration from sequence A at t i acquisition time to sequence B at t j acquisition time as follow:

M A,B ti,tj , i, j ∈ {0, 1}, A, B ∈ {T 1, Gado, F lair} 3.
Image segmentation: we segmented the structures of interest(i.e tumor, edema, ventricles, and the brain) on the images by combining manual segmentation using ITK-SNAP [START_REF] Yushkevich | User-guided 3D active contour segmentation of anatomical structures: significantly improved efficiency and reliability[END_REF] and automatic segmentation based on level set method. [START_REF] Ciofolo | Brain Segmentation with Competitive Level Sets and Fuzzy Control[END_REF] We noted by V A ti , i ∈ {0, 1}, A ∈ {T 1, Gado, F lair} the volume segmented from image sequence A at time acquisition t i .

4. Geometrical transformations of all volumes to the reference frame are noted as following :

V A,B ti,tj = V A ti * M A,B ti,ti * M B,B ti,tj
where V A,B ti,tj is the volume segmented from sequence A at acquisition time t j reformated to the sequence B at time t i , i, j ∈ {0, 1} and A, B ∈ {T 1, Gado, F lair}

Estimation of local deformations

The estimation of local deformations was divided into three parts :(1) introduction of a metric to evaluate the percentage of volume change and the nature of deformation (shrinking/dilation), (2) estimation of surface deformations of the structures components of the region of interest (ROI), and (3) estimation of the error introduced by the manual segmentation.

Shrinking/dilation metric

To estimate the volume change of different structures from acquisition time t 0 to time t 1 , we introduced a local volume change as . A positive (negative) value of ∆V represents, respectively, an increase (decrease) in the volume. The ∆V is a dilation/shrinking metric for results interpretation.

Distance maps visualization

To estimate the deformations of the surfaces of the segmented structures and to visualize these deformations according to the localization of the structures, we developed a tool to extract and visualize a signed Distance Maps (DM ) for the different surfaces. Theses surfaces correspond, respectively, to the segmented structures at time t 0 (reference structures S t0 ) and t 1 (target structures S t1 ). Without loss of generality, we assumed that 0 < DSC(S t0 , S t1 ) < 1, where DSC 14 (Dice Similarity Coefficient) measures the spatial overlap between the two segmentations (S t0 , S t1 ) and is defined as follow:

DSC(S t0 , S t1 ) = 2(S t0 ∩ S t1 ) (S t0 + S t1 )
The distance maps visualization was performed as following: Initialization:

1. Take two point sets S t0 and S t1 2. Initialize the distance maps DM (x, y) = 0 ∀x ∈ S t0 , y ∈ S t1 Iterations: 3. For each point x 0 ∈ S t0 , find the closest point x 1 ∈ S t1 4. Define the sign of the distance: sign(DM (x 0 , x 1 )) = (-)(i.e shrinking) if x 1 ∈ Interior(S t0 ), (+)(i.e dilation) else 5. Calculation of the distance: DM (x 0 , x 1 ) = sign(DM (x 0 , x 1 )) * x 1x 0 2 Visualization: 6. Display the distance maps for the surfaces S t0 and S t1 , with color coding proportional to the distance.

Error of manual segmentation and partial volumes effects

The expert * segmented the tumor twice for one case study at time t 0 and time t 1 . The variation of volume between both segmentations gives an approximation of the manual volume segmentation error. The manual segmentation is performed at voxel level. To evaluate the partial volume effects, the expert displayed one segmentation and delineated the corresponding contour. This 3 rd segmentation gives a relative volume change according to subvoxel ignorance.

Estimation of regional displacements

Landmarks displacement around the ROI

To quantify displacement of neighboring structures to the ROI, the expert identified eight anatomically homologous landmarks L k ti , i ∈ {0, 1}, k ∈ {0, ..8} using MINC tools † . We note L k 01 = (L k t0 , L k t1 ) the matching couple between landmark L k t0 at time t 0 and the corresponding landmark L k t1 at t 1 . We note by L k 01 2 = (L k t0 , L k t1 ) 2 the Euclidean distance between each k th landmark at time t 0 and the k th corresponding landmark at time t 1 .

Reference landmarks displacement vs registration

To validate the precision of landmarks identification according to the error of rigid-body registration used in the image processing WorkFlow, the expert selected eight additional landmarks that are not supposed to be deformed (e.g the nasion, the eyes). We note this set of reference landmarks F k ti , i ∈ {0, 1}, k ∈ {0, ..8}. The registration was be accepted if the mean of displacement of the different reference landmarks was less than 1mm.

Estimation of ROI components displacement: Center of mass displacement

We supposed that the displacement of the ROI components (Tumor, Edema) can be characterized by the movement of their center of mass. We defined center of mass for each volume V of ROI components (CM (V )) and we note (CM (V t0 ), CM (V t1 )) 2 the Euclidean distance between the volume V t0 at acquisition time t 0 and the volume V t1 at acquisition time t 1 .

Introduction of Surgical Planning (SP) Error

For each intervention, the SP contribution consists of two parts : (1) localization of the target, and (2) definition of a road map according to the anatomical references (cortical gyri or sulci) and constraints (functional areas). We defined the difference between SP based on the acquisition time t 0 (SP t0 ) and SP based on t 1 time (SP t1 ) by given couple of errors related to: (1) localization of targets (ROI) from t 0 to t 1 and (2) displacement of anatomical references (like cortical gyri or sulci) from t 0 to t 1 . If we classify the data according to the treatment, we can define the error as follow :

Error(SP t0 , SP t1 ) =
[ (target t 0 , target t 1 ) 2 ; (sulci t 0 , sulci t 1 ) 2 ] if resection, [ (target t 0 , target t 1 ) 2 ; M ean(landmarks displacement)] if biopsy.

Where : target t 0 , target t 1 are the targets (tumors for our three surgical cases) from acquisition time t 0 , t 1 respectively and sulci t 0 , sulci t 1 are the sulci from t 0 , t 1 .

RESULTS

Estimation of local deformations

Local volume changes

From Table 1, a significant increase in the tumor volumes for the first and the third patient contrasts with decrease of tumor volume for the second patient. However, the edema and the global ROI (edema + tumor) volumes decreased for all patients by at least 12.18%. The variation of volume of ventricles depends on tumor localization. We observe for the second patient a significant increase (+21.77% ).

Manual segmentation error

Table 2 shows, for the second patient, a non significant partial volume effects due to the manual segmentation by the expert. We observed for the two segmentations of the tumor from time t 1 a variation: ∆S T umor = -0.407%.

Comparing to the volume changes of the tumor for the second patient (Table 1), we observed:

∆S T umor = -0.407% -28.60% = ∆V T umor

Distance maps visualization

From Figure 2, we distinguished the shrinking/dilation of different surfaces of structures from time acquisition t 0 to time t 1 . The visualization of these distance maps shown location of the shrinking or dilation and illustrates significant deformations of surfaces and confirmed the results of volume changes. The visualization of brain distance map from t 0 to t 1 merged with ROI localization shows a shift of the brain according to the shrinking/dilation of the structures component of the ROI. Volumes for tumor segmentation at t 0 and t 1 S1(mm From top to bottom, the three rows correspond to the 1 st , 3 rd , and 2 nd patient respectively. The last row represents for the 2 nd patient the merging between the tumor localization and the brain distance map.

Estimation of regional displacement

Table 3 shows a significant displacement of landmarks around the ROI of the three patients. For the 2 nd patient, the mean of displacement was about 3.90 ± 0.92 mm with maximal displacement of 5.45 mm. In second part, the table shows that the rigid registration applied for our IPWF was accepted. For the three patients, the means were: 0.54±0.31 mm, 0.63±0.30 mm, and 0.83±0.22 mm respectively. The mean of distances between matching reference landmarks was inferior to 1 mm.

Estimation of ROI component displacement

From (Table 4) we observed a significant displacement of mass center of the ROI (CM(ROI)); for the 3 rd patient, the CM(edema+tumor) moved by 6.92 mm and for the 1 st patient, the center of mass of tumor moved by 4.05 mm. These results coupled with landmarks displacement demonstrate first sizable difference between SP based t 0 vs SP based t 1 .

Surgical Planning (SP) Error

For the second patient, the surgeon defined an optimal sulci, to reach the target, for the planning SP t0 and the corresponding sulci for the planning SP t1 . The left part of the (Figure 3) shows the difference between the sulci extraction based on time t 0 (yellow color) and the extraction based on time t 1 (blue color). The right part of the (Figure 3) shows the optimal reference sulci to access to the tumor based on SP t0 (yellow) and the corresponding one according to SP t1 (blue color). The Euclidean distance between the center of mass of the sulci from t 0 to t 1 was (sulci t 0 , sulci t 1 ) 2 = 7.95 mm and from (Table 3), we note (tumor t 0 , tumor t 1 ) 2 = 2.70 mm. It gives a surgical planning error for the second patient as:

Error(SP t0 , SP t1 ) = [2.70; 7.95]mm

We resume in Table 5 the corresponding errors due to planning based on t 0 vs planning based on t 1 for the three patients.

CONCLUSION AND DISCUSSION

Different studies have focused on compensating the brain shift phenomenon and to correct deformation in intraoperative stage. [START_REF] Roberts | Intraoperative brain shift and deformation: a quantitative analysis of cortical displacement in 28 cases[END_REF][START_REF] Dumpuri | An atlas-based method to compensate for brain shift: Preliminary results[END_REF] Since temporal aspect is a key factor during neurosurgery, to extract an information about deformation can increase the accuracy of the MIGNS systems. Different parameters related to the phenomenon increase the complexity of the event. Drug responses, gravity, edema and surgical manipulation complicate the understanding of the deformation. [START_REF] Miga | Model-updated image guidance: initial clinical experiences with gravity-induced brain deformation[END_REF][START_REF] Dumpuri | Automated brain shift correction using a pre-computed deformation atlas[END_REF] The brain shift phenomenon is certainly related to the disturbance of the volume equilibrium and a studies [START_REF] Roberts | Intraoperative brain shift and deformation: a quantitative analysis of cortical displacement in 28 cases[END_REF][START_REF] Hill | Measurement of introperative brain surface deformation under a craniotomy[END_REF] have shown a cortical brain surface shift of 10 mm on average, during surgery and along direction of deformations, after skull opening.

In this study, we addressed another aspect of anatomical deformation related to brain shift in the preoperative phase that we call preoperative brain shift. This preoperative brain shift was restricted to the anatomical brain deformations in closed-cranium. We hypothesized, for our data set of patients (three surgical cases), that the presence of edema and tumor creates a mass effect and that after corticosteroids administration during stage [t 0 , t 1 ], both the shrinking of edema and growth of tumor induce lack of tissue and provide displacement of neighboring structures of the ROI. The three surgical cases illustrate the existence of the preoperative brain shift and initialize preliminary measures to understand the shift in closed cranium under controlled parameters (tumor localization, edema volume, tumor volume, and corticosteroids administration) without surgical manipulation. The signed distance maps introduced in this study can be used for understanding and modeling preoperative deformation. According to brain surface deformations, we will be able to localize an optimal craniotomy in order to minimize the shift due to the decrease of intracranial pressure after skull opening. The local volume changes of the components of the ROI can be related to the localization of these components; the significant increase of the volume of ventricles for the second patient (Table 1) can be explained by the fact that the tumor and the edema are inside the ventricles.

On the methodological side, our study suffers from three principal points: (1) rigid body registration, (2) manual segmentation, and (3) intra-expert variability in landmarks identification. For further work, we will introduce elastic body registration to extract the deformation field and to identify the features for direction of deformation. For manual segmentation by the expert, we observed that we can ignore the partial volume effects comparing to the significant volume changes, but we have to introduce the error related to different experts. For homologous landmarks identification, we have to estimate the displacement errors associated with placement of the landmarks and errors of locating the fiducial points due to the intra-expert variability. Another further step of our work will be to propose a realistic model of tumor growth, that is in general difficult, [START_REF] Clatz | Realistic Simulation of the 3D Growth of Brain Tumors in MR Images Coupling Diffusion with Mass Effect[END_REF] and combine it with model of edema shrinking in order to build a first model for preoperative brain shift.

Further analysis, with more clinical cases, will certainly emphasize a new surgical workflow that determine a time t 1 ∈ [t 0 , t S [ (t S time of intervention), respecting clinical temporal constraints and processing time of images (preprocessing, fusion and segmentation), after which new preoperative imaging must be performed in order to initialize properly the MIGNS systems and recover the error of misregistration. In fine, this would be a proper initialization of the intraoperative brain deformation and help approach the reality in the operating room.

Figure 1 .

 1 Figure 1. Image processing workflow for volumes segmentation and registration of all volumes in the same reference.

Vt 1

 1 Vt 0 and the relative volume change ∆V = Vt 1 -Vt 0 Vt 0

Figure 2 .

 2 Figure 2. Signed distance maps, based on closest points algorithm, from reference time t0 to targeted time t1. From left to right, distance maps for tumor, ROI(Edema+Tumor), and ventricles.From top to bottom, the three rows correspond to the 1 st , 3 rd , and 2 nd patient respectively. The last row represents for the 2 nd patient the merging between the tumor localization and the brain distance map.

Figure 3 .

 3 Figure 3. The extraction of sulci based on time t0 with optimal reference sulci (yellow color) and the extraction based on time t1 (blue color) with optimal reference sulci to access to the tumor.

Table 1 .

 1 Table representing volume changes ∆V for the three patients

	Structures / patient	Volume changes ∆V = V1-V0 V0 % 1 st patient 2 nd patient 3 rd patient
	Brain	-4.92%	+1.64%	+2.01%
	Tumor	+60%	-28.60%	+16.58%
	Edema	-32.92%	-31.00%	-15.34%
	ROI=Edema+Tumor	-13.62%	-30.74%	-12.18%
	Ventricles	-0.94%	+21.77%	-13.87%

Table 2 .

 2 Error of manual volume segmentations for tumor (2 nd patient)

Table 3 .

 3 Euclidean distance . 2 between each landmark L k t 0 at time t0 and the corresponding landmark L k t 1 at time t1 for the three patients.

		Displacement of Landmarks around the ROI
	Patient/ . 2 (mm) L 1 01 1 st patient 2.14 0.84 1.17 2.36 1.87 1.09 0.89 2.77 1.64 0.69 L 2 01 L 3 01 L 4 01 L 5 01 L 6 01 L 7 01 L 8 µ σ 01
	2 nd patient	2.08 3.07 4.17 4.13 4.36 4.25 3.68 5.45 3.90 0.92
	3 rd patient	2.16 1.47 1.78 2.47 1.73 2.65 1.87 2.89 2.13 0.46
		Displacement of reference landmarks
	Patient/ . 2 (mm) F 1 01 1 st patient 0.17 0.25 1.15 0.63 0.43 0.93 0.44 0.34 0.54 0.31 F 2 01 F 3 01 F 4 01 F 5 01 F 6 01 F 7 01 F 8 µ σ 01
	2 nd patient	0.37 0.30 0.43 0.72 1.27 0.89 0.49 0.61 0.63 0.30
	3 rd patient	0.93 1.19 0.66 0.78 0.87 0.40 0.80 1.01 0.83 0.22

Table 4 .

 4 The displacement of Center of Mass (CM) for the tumor and the global ROI (Edema+Tumor)

		CM displacement of The ROI components
	Patient/ . 2 (mm) CM(Tumor) 1 st patient 4.052	CM(ROI) 1.713
	2 nd patient	2.706	1.45
	3 rd patient	1.971	6.920

Table 5 .

 5 Surgical planning error for the three patients

	Surgical planning error
	1 st patient 2 nd patient 3 rd patient
	Error(SP t0 , SP t1 )mm [4.05; 1.64] [2.70; 7.95] [1.97; 2.13]

ACKNOWLEDGMENTS

We wish to express our thanks to Pr. G. Brassier and all people of the department of neurosurgery (C. Haegelen, S. Diabira) and radiology (J-C. Ferré, B. Carsin) of the Pontchaillou Hospital. We thank N. Wiest-Daesslé and S. Prima for review. Our research is supported by FRSQ/INSERM (S. Duchesne) and a french research ministry grant (O. El Ganaoui).


Ahmed Zeriahi
email: ahmed.zeriahi@math.univ-toulouse.fr

REMARKS ON THE MODULUS OF CONTINUITY OF SUBHARMONIC FUNCTIONS

Keywords: . 2010 Mathematics Subject Classification. 31B05, 31C05, 32U05 Subharmonic functions, Poisson formula, Modulus of continuity, quasi-plurisubharmonic functions

We introduce different classical characteristics used to regularize a subharmonic function and compare them.

As an application we give a complete proof of a useful characterization of the modulus of continuity of such functions in terms of these characteristics under a technical condition. This result is extended to quasi-plurisubharmonic functions on a compact Hermitian manifold.

Introduction

Given a subharmonic function u on a domain Ω ⊂ R n , we introduce various continuous approximating functions of u and use them to define various "partial moduli of continuity" associated to u. These moduli have been used in many papers to measure in different ways the continuity of solutions to complex Monge-Ampère equations on bounded domains in C n as well as on compact complex manifolds (see [START_REF] Guedj | Hölder continuous solutions to Monge-Ampère equations[END_REF], [START_REF] Demailly | Hölder Continuous Solutions to Monge-Ampère Equations[END_REF], [N18], [START_REF] Ko Lodziej | A remark on the continuous subsolution problem for the complex Monge-Ampère equation[END_REF], [START_REF] Benali | The Hölder continuous subsolution theorem[END_REF]) .

The goal of this note is to clarify the relations between these moduli and establish estimates on the (full) modulus of continuity of a subharmonic function in terms of these partial moduli of continuity.

Let Ω ⊂ R n be a domain and u : Ω -→ R ∪ {-∞} be a subharmonic function on Ω.

We fix δ 0 > 0 so that Ω δ 0 := {x ∈ Ω ; dist(x, ∂Ω) > δ 0 } = ∅. For 0 < δ < δ 0 , we can define the δ-mean value function associated to u as follows We can also consider the δ-max regularization of u defined as follows:

(1.2) M δ u(x) := max y∈B u(x + δy), for x ∈ Ω δ and 0 < δ < δ 0 .

We have obviously Λ δ u(x) ≤ M δ u() for any x ∈ Ω δ . Moreover these functions are continuous and subharmonic on Ω δ and decrease to u pointwise on Ω.

We consider a modulus of continuity κ : R + → R + which satisfies the following condition:

(1.3) ∃A > 0, lim sup

t→0 + κ(At) Aκ(t) < 1 2n •
The main result of this note is the following.

Main Theorem. Let u : Ω -→ R be a bounded subharmonic function on Ω and κ a modulus of continuity satisfying the condition (3.1). Assume that there exists a constant C 0 > 0 and 0 < δ 1 < δ 0 such that for 0 < δ < δ 1 ,

(1.4) Λ δ u(x) -u(x) ≤ C 0 κ(δ), for x ∈ Ω δ .
Then there exists constants B > 1, ε 0 > 0 with Bε 0 < δ 1 and a constant C > 0 such that for 0 < δ < ε 0 and x ∈ Ω Bδ , we have

(1.5) M δ u(x) -u(x) ≤ C κ(δ).
In particular u is κ-continuous on any compact subset E Ω.

Observe that the condition of the theorem is satisfied for any harmonic function, regardless on its regularity at the boundary. Hence we cannot expect to conclude anything about the full modulus of continuity of u on the whole domain Ω.

However if we know about the behaviour of u near the boundary we can get a better control on the modulus of continuity of u on Ω.

Corollary. Let u : Ω -→ R be a bounded subharmonic function on Ω and κ a modulus of continuity satisfying the condition (1.3).

Assume that u satisfies the condition (1.4) and extends as a κ-continuous function near the boundary. Then u is κ-continuous on Ω i.e.

(1.6) |u(x) -u(y)| ≤ L κ(|x -y|),
for any x, y ∈ Ω, where L > 0 is a uniform constant.

Let us recall the definition of κ-continuity near the boundary. Set for

δ > 0, κu (δ) := sup{|u(x) -u(y)|; y ∈ ∂Ω, x ∈ Ω ∩ B(y, δ)}•
It is easy to see that if u is continuous on ∂Ω, then lim δ→0 κu (δ) = 0. Then we say that u : Ω -→ R is κ-continuous near the boundary if there exists a constant C > 0 such that κu (δ) ≤ Cκ(δ) for δ > 0 small enough.

Characteristics associated to a subharmonic functions

Let Ω ⊂ R n be a domain. We denote by SH(Ω) ⊂ L 1 loc (Ω) the set of subharmonic functions on Ω. We will will first introduce differents characteristics associated to u and compare them. Then we prove some average estimates on them.

2.1. Basic definitions. Let u ∈ SH(Ω). We will associate to u the following characteristics. Set Ω δ := {x ∈ Ω; dist(x, ∂Ω) > δ} for 0 < δ ≤ δ 0 , where δ 0 > 0 is fixed so that Ω δ 0 = ∅.

Fix 0 < δ < δ 0 and x ∈ Ω δ . Define the max characteristic of u as follows

(2.1) M δ u(x) := max B(x,δ) u = max |ξ|=1 u(x + δξ),
where B(x, r) := {y ∈ R n ; |y -x| ≤ r} is the closed euclidean ball of center x and radius r > 0.

We define the mean value volume characteristic of u

(2.2) Λ δ u(x) := 1 τ n r n B(x,δ) u(y)dλ n (y) = 1 τ n B u(x + δy)dλ n (y),
where dλ n is the Lebesgue measure on R n and τ n := λ n (B) is the volume of the unit ball B ⊂ R n . We define the mean value area characteristic of u

(2.3) A δ u(x) := 1 σ n-1 S u(x + δξ)dσ(ξ),
where dσ is the area measure of the unit sphere S = ∂B and σ n-1 is the area of the unit sphere S ⊂ R n . We consider more general mean value characteristics associated to the subharmonic function u.

Let ρ be a radial bounded Borel function with compact support in the unit ball B ⊂ R N such that B ρ(x)dλ n (x) = 1. This means in spherical coordinates that

(2.4) σ n-1 . 1 0 ρ(r)r n-1 dr = 1.
For δ > 0 and x ∈ R n , we set ρ δ (x) := δ -n ρ(x/δ). Then it's well known that ρ δ → 0 , in the sense of distributions on R n as δ → 0, where 0 is the unit mass Dirac distribution at the origin.

For x ∈ Ω δ , the smooth mean value characteristic of u is defined by (2.5)

R δ u(x) = u ρ δ (x) := B(x,δ) u(ξ)ρ δ (x-ξ)dλ n (ξ) = B u(x+δξ)ρ(ξ)dλ n (ξ).
Observe that if ρ = 1 τn 1 B then R δ u(x) = Λu δ (x). It is well known that all these functions are non decreasing in δ and convex in the variable t = k n (r), where k 2 (r) = log r and k n (r) := -r 2-n when n ≥ 3 (see [START_REF] Armitage | Classical Potential Theory[END_REF]).

2.2. Comparison of characteristics. We want to compare all these characteristics for a subharmonic function.

Lemma 2.1. Let Ω ⊂ R n be an open set and u ∈ SH(Ω). Then for any 0 < δ < δ 0 and x ∈ Ω δ , we have

b n A δ/2 u(x) -u(x) ≤ R δ u(x) -u(x) ≤ A δ u(x) -u(x), In particular, b n A δ/2 u(x) -u(x) ≤ Λu δ (x)-u(x) ≤ A δ u(x)-u(x) ≤ M δ u(x)-u(x),
where b n := 1 1/2 ρ(r)r n-1 dr < 1/σ n-1 . Proof. Fix 0 < δ < δ 0 and x ∈ Ω δ . Using spherical coordinates, we see that

R δ u(x) = 1 0 ρ(r)r n-1 S u(x + rδξ)dσ(ξ) dr = σ n-1 1 0 ρ(r)r n-1 A rδ u(x)dr.
Therefore by the equation (2.4), we deduce that for any fixed

x ∈ Ω δ (2.6) R δ u(x) -u(x) = σ n-1 1 0 ρ(r)r n-1 (A rδ (x) -u(x)) dr.
Since the function ]0, δ] s → A s u(x) -u(x) in non-negative and nondecreasing, it follows that

(2.7) b n (A δ/2 u(x) -u(x)) ≤ R δ u(x) -u(x) ≤ A δ u(x) -u(x),
where b n := 1 1/2 ρ(r)r N -1 dr < 1. Now applying the formula (2.6) with ρ = 1 τn 1 B , we obtain

Λ δ/2 u(x) -u(x) = n 1 0 r n-1 (A rδ/2 u(x) -u(x))dr, since σ n-1 = nτ n .
Applying the inequality (2.7) we obtain for 0 < δ < δ 0 and x ∈ Ω δ

Λ δ/2 u(x) -u(x) ≤ A δ/2 u(x) -u(x) ≤ b -1 n (R δ u(x) -u(x)).
Now we want to compare the supnorm and the mean value of a subharmonic function on balls.

Lemma 2.2. There exists δ 0 > 0 small enough and a constant a n > 0 such that for any 0 < δ < δ 0 , 0 < θ < 1 and x ∈ Ω δ , we have

M θδ u(x) -u(x) ≤ c n (A δ u(x) -u(x)) + c n 2 n θ (1 -θ) n-1 S (u(x + δy) -u(x)) + dσ(y)• (2.8)
In particular

M δ/2 u(x) -u(x) ≤ c n (A δ u(x) -u(x)) + 4 n c n S (u(x + δy) -u(x)) + dσ(y).
(2.9)

Proof. Assume first that u is subharmonic in a neighbourhood of the closed ball B. It follows from Poisson-Jensen formula for the unit ball B (see [START_REF] Armitage | Classical Potential Theory[END_REF]) that (2.10) u(x) ≤ S P (x, y)u(y)dσ(y), x ∈ B, where

P (x, y) := c n 1 -|x| 2 |x -y| n , (x, y) ∈ B × ∂B, is the Poisson kernel of the unit ball B ⊂ R n .
Since S P (x, y)dσ(y) = 1, it follows from (2.10) that for any x ∈ B,

(2.11) u(x) -u(0) ≤ S P (x, y)(u(y) -u(0))dσ(y) = I + (u) + I -(u),
where (2.12)

I + (u) := {y∈S;u(y)≥u(0)} P (x, y)(u(y) -u(0)) dσ(y),

and

(2.13)

I -(u) := {y∈S;u(y)≤u(0)} P (x, y)(u(y) -u(0)) dσ(y).
Since for 0 < r < 1 and |x| = r, we have

(2.14) c n 1 -r (1 + r) n-1 ≤ P (x, y) ≤ c n 1 + r (1 -r) n-1 ,
it follows from (2.12) and (2.14) that for |x| = r < 1, (2.15)

I + (u) ≤ c n 1 + r (1 -r) n-1 {y∈S;u(y)≥u(0)} (u(y) -u(0)) dσ(y).
Moreover, it follows from (2.13) and (2.14) that for |x| = r < 1, (2.16)

I -(u) ≤ c n 1 -r (1 + r) n-1 {y∈S;u(y)≤u(0)} (u(y) -u(0)) dσ(y).
Therefore from (2.15), (2.16) and (2.11), we deduce that

u(x) -u(0) ≤ c n 1 + r (1 -r) n-1 {y∈S;u(y)≥u(0)} (u(y) -u(0)) dσ(y) + c n 1 -r (1 + r) n-1 {y∈S;u(y)≤u(0)} (u(y) -u(0)) dσ(y).
Observe that for r < 1,

1 + r (1 -r) n-1 ≤ 1 + 2 n r (1 -r) n-1 , and 1 -r (1 + r) n-1 ≥ 1 -nr. This implies that max |x|=r u(x) -u(0) ≤ c n S (u(y) -u(0)) dσ(y) + c n 2 n r (1 -r) n-1 S (u(y) -u(0)) + dσ(y) -nc n r S (u(y) -u(0)) -dσ(y).
Now in the general case fix δ 0 > 0 small enough so that Ω 2δ 0 = ∅. We fix x ∈ Ω 2δ so that B(x, 2δ) ⊂ Ω and apply the previous inequality to the function y -→ u(x + δy) which is subharmonic in a neighbourhood of the unit ball B. We then obtain for 0 < r < 1 max

B(x,rδ) u -u(x) ≤ c n S (u(x + δy) -u(x)) dσ(y) + c n 2 n r (1 -r) n-1 S (u(x + δy) -u(x)) + dσ(y) -nc n r S (u(x + δy) -u(x)) -dσ(y).
Since the last term is non positive, the inequality of the lemma follows.

We can easily deduce the following result.

Corollary 2.3. Let u be a bounded subharmonic function on a bounded domain Ω ⊂ R n . Assume that there exists α ∈]0, 1] and κ 1 > 0 such that for any 0 < δ < 2δ 0 and x ∈ Ω δ , we have

R δ u(x) ≤ u(x) + κ 1 δ α .
Then there exists κ 2 > 0 and 0 < δ 1 << 1 such that for any δ < δ 1 and x ∈ Ω δ , we have

M δ u(x) ≤ u(x) + κ 2 δ α/(1+α) .
Proof. Apply Lemma 2.2 with θ = δ α . Then for δ < δ 1 < 2 -1/α and x ∈ Ω δ 1+α , we have

M δ 1+α u(x) -u(x) ≤ c n κ 1 δ α + L n M δ α .
where M := oscΩu is the oscillation of u on Ω and L n > 0 is a uniform constant. Relpacing δ by δ 1/(1+α) we obtain that u is Hölder continuous with exponent α/(1 + α) and κ 2 := a n κ 1 + nM L n .

2.3. Average estimates. We first recall a well know result which is important in applications. This result is shown in [?], beut we will recall the proof here for the convenince of the reader.

Lemma 2.4. Let u be subharmonic function on a bounded domaine Ω ⊂ R n . Then there exists a uniform constant e n > 0 such that for 0 < δ < δ 0 ,

Ω δ (A δ u(x) -u(x)) dy ≤ e n ∆u Ω δ δ 2 .
Proof. Let µ := (1/2π)∆u be the Riesz measure of u on Ω. It follows from Poisson-Jensen formula that for x ∈ Ω with u(x) > -∞, we have

A δ u(x) -u(x) = δ 0 t 1-n (µ(B(x, t))dt.
Then integrating on x over Ω and applying Fubini's Theorem we obtain

A δ u(x) -u(x)dλ n (x) ≤ δ 0 t 1-n Ω δ dµ(ζ) (B(ζ,t) dλ n (z)dt ≤ τ n µ(Ω δ ) δ 0 t dt = τ n µ(Ω δ )δ 2 /2,
where τ n = λ n (B) is the volume of the unit ball in R n . This proves the Lemma.

Corollary 2.5. Assume that u is a subharmonic function on a bounded domain Ω ⊂ R n . Then for any 0 < δ < δ 0 /2, we have

Ω 2δ (M δ u -u)dλ n ≤ p n Ω 2δ (A 2δ u -u)dλ n + q n δ ∇u L 1 (Ω δ ) ,
where p n , q n > 0 are uniform constants.

In particular

Ω 2δ (M δ u -u)dλ n ≤ p n ∆u Ω 2δ δ 2 + q n ∇u L 1 (Ω δ ) δ,
where p n , q n > 0 are uniform constants.

Proof. We apply Lemma 2.2 and integrate over Ω 2δ . Then for 0 < δ < 2δ 0 ,

Ω 2δ (M δ u(x) -u(x))dλ(x) ≤ c n b n Ω 2δ (A 2δ u(x) -u(x)) + c n 2 2n-2 S Ω 2δ |u(x + δy) -u(x)|dλ(x) dσ(y),
where p n = c n b n and q n := c n 2 2n-2 . We claim that for fixed x ∈ Ω δ and 0 < δ < δ 0 /2, we have

(2.17)

Ω 2δ |u(x + δy) -u(x)|dλ(x) ≤ δ ∇u L 1 ( Ωδ ) .
Indeed assume first that u is smooth. Now observe that for |y| = 1, we have

u(x + δy) -u(x) = δ 0 Du(x + ty) • ydt.
Then using Fubini's Theorem, we obtain Integration over Ω 2δ leads to the inequality (2.17). For a non smooth function u, we can approximate u by a decreasing sequence (u j ) of smooth subharmonic functions on a neighbourhood V of Ωδ so that Du j → Du in L 1 (V ) and almost everywhere on V . Applying the inequality (2.17) to the u j 's and passing to the limit, we obtain by Fatou's lemma the inequality (2.17) for u.

More general moduli of continuity

We want to show that the sup-regualrization and the mean value regularisation have the same behaviour for a large class of moduli of continuity. Namely we will prove an important result which confirms a lemma stated in [START_REF] Guedj | Hölder continuous solutions to Monge-Ampère equations[END_REF] and used in the litterature for a Hölder modulus of continuity. Chinh H. Lu discovered recently a gap in the proof of [START_REF] Guedj | Hölder continuous solutions to Monge-Ampère equations[END_REF] which was fixed in [START_REF] Lu | Stability and Hölder regularity of solutions to complex Monge-Ampère equations on compact hermitian manifolds[END_REF] in the case of a compact hermitian manifold (without boundary) following the same scheme.

We will follow the same scheme as in [START_REF] Guedj | Hölder continuous solutions to Monge-Ampère equations[END_REF] and use a new idea of [START_REF] Lu | Stability and Hölder regularity of solutions to complex Monge-Ampère equations on compact hermitian manifolds[END_REF] to prove a more general result.

A new characterization. Let us first give some definitions.

Let κ : [0, l] -→ R + be a modulus of continuity i.e. a continuous increasing subadditive function such that κ(0) = 0.

We will consider the following growth condition on κ.

(3.1) ∃A > 0, lim sup

t→0 + 2nκ(At Aκ(t) < 1.
Observe that this condition holds for a logarithmic Hölder modulus of continuity defined by κ α,β (t) = t α (-log t) β , 0 < t < t 0 < 1 with 0 ≤ α < 1 and β ∈ R, with β < 0 if α = 0 and t 0 > 0 is chosen so small that κ α,β is concave on [0, t 0 ]. However it's not satisfied by the modulus of continuity κ 1,β with β ≤ 0.

We need another definition.

Definition 3.1. We say that a function u : Ω -→ R is κ-continuous near the boundary ∂Ω if there exits ε 0 > 0 small enough such that for any ζ ∈ ∂Ω and z ∈ Ω with |z -ζ| ≤ ε 0 , we have

(3.2) |u(z) -u(ζ)| ≤ κ(|z -ζ|).
Observe that this condition implies the continuity of u on ∂Ω and it is satisfied if there exists two functions v, w κ-continuous near the boundary such that v ≤ u ≤ w near the boundary and v = u = w on ∂Ω.

We need to introduce one more characteristic associated to u. For 0 < δ < δ 0 , we set

O δ u(x) := osc B(x,δ) u = max{|u(y 1 ) -u(y 2 )|; y 1 , y 2 ∈ B(x, δ)}•
Now we can sate the main result of this note. Theorem 3.2. Let κ be a modulus of continuity satisfying the condition (3.1) and u : Ω -→ R be a bounded function which is subharmonic on Ω and κ-continuous near the boundary ∂Ω.

Then the following properties are equivalent :

(1) there exists a constant L 1 > 0 and 0 < δ 1 < δ 0 such that for

0 < δ < δ 1 , M δ u(x) -u(x) ≤ L 1 κ(δ), for x ∈ Ω δ , (2
) there exists constants L 2 > 0 and 0 < δ 2 < δ 0 such that Bδ 2 < δ 0 and for 0 < δ < δ 2 ,

R δ u(x) -u(x) ≤ L 2 κ(δ), for x ∈ Ω δ ,
(3) there exists a constants B > 1, L 3 > 0 and 0 < δ 3 < δ 0 such that for 0 < δ < δ 3 , O δ u(x) ≤ L 3 κ(δ), for x ∈ Ω Bδ , (4) the function u is κ-continuous on Ω i.e. there exists a constant L 3 > 0 and 0 < δ 3 < δ 0 such that for any x ∈ Ω and y ∈ Ω with |x -y| ≤ δ 3 , we have

|u(x) -u(y)| ≤ L 3 κ(|x -y|).
Observe that the condition (2) is always satisfied for any harmonic function u on Ω, regardless of its behaviour at the boundary, while the condition (4) implies that the boundary values of u is κ-continuous on ∂Ω. Therefore (2) and ( 4) are not equvalent without any condition on the behaviour of u at the boundary.

The main step in the proof of our theorem is the following lemma whose proof is inspired from [START_REF] Guedj | Hölder continuous solutions to Monge-Ampère equations[END_REF] and [START_REF] Lu | Stability and Hölder regularity of solutions to complex Monge-Ampère equations on compact hermitian manifolds[END_REF].

Lemma 3.3. Let u : Ω -→ R be a bounded subharmonic function on Ω and κ a modulus of continuity satisfying the condition (3.1) . Assume that there exists a constant C 0 > 0 and 0 < δ 1 < δ 0 such that for 0 < δ < δ 1 ,

(3.3) R δ u(x) -u(x) ≤ C 0 κ(δ), for x ∈ Ω δ .
Then there exists constants B > 1, ε 0 > 0 with Bε 0 < δ 1 and a constant C 3 > 0 such that for 0 < δ < ε 0 and x ∈ Ω Bδ , we have

(3.4) O δ u(x) ≤ C 3 κ(δ).
In particular u is κ-continuous on any compact set E Ω.

Proof. We claim that the condition (3.1) implies that we can choose A > 2 large enough and δ 2 > 0 small enough such that (2A + 1)δ 2 < δ 1 and

(3.5) θ := 2n sup

0<t≤δ 2 κ((A + 1)t) Aκ(t) < 1.
Indeed by (3.1) we see that there exists A > 0, ν < 1/2n and 0 < δ 2 < δ 0 small enough such that sup

0<s≤δ 2 κ(A s) A κ(s) < ν < 1/2n.
Fix an integer N > 1 and apply this inequality for s = N t. Then by subadditivity of κ, we have for 0

< t < δ 2 := δ 2 /N , κ((N A t) N A κ(t) ≤ κ(N A t) A κ(N t) = κ(A s) A κ(s) < ν. Now choose N > 1 so large that N A > 3 and set A := N A -1 > 2.
Then the previous inequality implies that for 0 < t < δ 2 , κ(((A + 1)t)

Aκ(t) = κ((N A t) N A κ(t) N A A < ν N A A .
Since N A A = N A N A -1 → 1 as N → +∞ and ν < 1/2n, we can find N > 2 large enough so that ν N A A < 1/2n, which implies the inequality (3.5) and proves the claim. Now choose δ 2 > 0 so small that (2A + 1)δ 2 < δ 1 and fix 0 < δ < δ 2 . The first step of the proof follows the scheme given in [START_REF] Guedj | Hölder continuous solutions to Monge-Ampère equations[END_REF]. Observe that the condition (3.3) implies that u is continuous on Ω.

Let x 0 ∈ Ω (2A+1)δ . Then B(x 0 , δ ⊂ Ω and by continuity, there exists ξ 0 , y 0 ∈ B(x 0 , δ) such that

(3.6) O δ u(x 0 ) = u(y 0 ) -u(ξ 0 ).
Since dist(x 0 , ∂Ω) > (2A + 1)δ, we have B(ξ 0 , 2Aδ) ⊂ Ω. Set a := A -2 > 0 and r := Aδ = (a + 2)δ. Then by Lemma 2.1, it follows that

(3.7) R 2r u(ξ 0 ) -u(ξ 0 ) ≥ b n (Λ r u(ξ 0 ) -u(ξ 0 )) , where b n := σ n-1 1 1/2 s n-1 ρ(s)ds. Since B(y 0 , aδ) ⊂ B(ξ 0 , r) ⊂ Ω, we can write, 1 τ n r n B(ξ 0 ,r) u(y)dλ n (y) = 1 τ n r n B(y 0 ,aδ) u(y)dλ n (y) + 1 τ n r n B(ξ 0 ,r)\B(y 0 ,aδ) u(y)dλ n (y). (3.8)
By subharmonicity of u, we have 1

τ n r n B(y 0 ,aδ) u(y)dλ n (y) ≥ a n (a + 2) n u(y 0 ).
On the other hand, we have

1 τ n r n B(ξ 0 ,r)\B(y 0 ,aδ) u(y)dλ n (y) ≥ 1 - a n (a + 2) n (u(y 0 ) -O r+δ u(x 0 )) .
Therefore adding the last two inequalties, we obtain

Λ r u(ξ 0 ) -u(ξ 0 ) ≥ a n (a + 2) n u(y 0 ) -u(ξ 0 ) + 1 - a n (a + 2) n (u(y 0 ) -O r+δ u(x 0 )) = O δ u(x 0 ) -1 - a n (a + 2) n O r+δ u(x 0 ).
From (2.11) and the previous estimate, we deduce that

R 2r u(ξ 0 ) -u(ξ 0 ) ≥ b n (O δ u(x 0 ) -a n O r+δ u(x 0 )) ,
where a n := (1 -a n (a+2) n . The second step of the proof will use an idea of [START_REF] Lu | Stability and Hölder regularity of solutions to complex Monge-Ampère equations on compact hermitian manifolds[END_REF]. Set f (t) := O t u(x 0 ) and g(t) := R 2t u(ξ 0 ) -u(ξ 0 ) for 0 < t < δ 2 and observe that a n ≤ 2n A , with A := a + 2. Then for 0 < δ < δ 2 , we have

f (δ) ≤ (2n/A)f ((A + 1)δ) + (1/b n )g(Aδ).
Now since ξ 0 ∈ Ω 2Aδ and 2Aδ < δ 1 , it follows from (3.3) that g(Aδ) ≤ C 0 κ(2Ar). Then by subadditivity of κ, we have for 0

< δ < δ 2 f (δ) ≤ (2n/A)f ((A + 1)δ) + C 1 κ(δ),
where

C 1 := (2A + 1)C 0 /b n .
Consider the quotient function h(t) := f (t)/κ(t). From the previous estimate we deduce that for 0 < δ < δ 2 ,

h(δ) ≤ 2n κ((A + 1)δ) Aκ(δ) h(Bδ) + C 1 ,
where B := A + 1. Now as A > 2 is choosen so that the condition (3.1) is satisfied, there exits 0 < δ 3 < δ 2 so that for 0 < δ < δ 3 , we have 2n

=: θ < 1. Then for 0 < δ < δ 3 we have

h(δ) ≤ θh(Bδ) + C 1 ,
Iterating this inequality we see that for any 0 < δ < δ 3 and any k ∈ N, we have

(3.9) h(δB -k ) ≤ θ k h(δ) + C 2 ,
where C 2 := C 1 1 1-θ . We claim that this inequality implies that h(t) is bounded near 0. Indeed choose ε 0 > 0 such that B 2 ε 0 < δ 3 and set

M 0 := max{h(t) ; ε 0 ≤ t ≤ B 2 ε 0 }.
Fix 0 < t < ε 0 and choose an integer k so that tB k ∈ [ε 0 , B 2 ε 0 ]. Then applying the inequality (3.9) with δ = tB k we obtain

h(t) ≤ θ k h(tB k ) + C 2 ≤ M 0 + C 2 =: C 3 ,
which proves our claim.

Let E Ω be a compact set and r 0 := dist(E, ∂Ω) > 0. Then for 0 < δ < r 0 /(2A + 1) , E ⊂ Ω (2A+1)δ . Then we can apply the previous estimate and get for x ∈ E and y ∈ E with |x -y| ≤ δ, u(x) -u(y) ≤ C 3 κ(δ).

We do not know if the condition (3.3) implies that the condition (3.4) holds for any x ∈ Ω δ , nor if it implies that u is κ-continuous on Ω.

3.2. Proof of the main theorem. We are now ready to prove the main theorem stated in the introduction.

Proof. It follows from from the assuptions od the main theorem that u is continuous on Ω. Therefore we can find x 0 ∈ Ω, ξ 0 ∈ Ω such that |x 0 -ξ 0 | ≤ δ and (3.10)

κ u (δ) = sup ξ,x∈ Ω(u(ξ) -u(x)) = u(ξ 0 ) -u(x 0 ).
Take 0 < δ 3 small enough so that (B + 1)δ 3 < δ 2 and fix 0 < δ < δ 3 . Then there are two cases to be considered for the point x 0 .

1) If x 0 ∈ Ω and dist(x 0 , ∂Ω) > Bδ, then x 0 ∈ Ω Bδ and ξ 0 ∈ B(x 0 , δ) and then by the inequality (3.4) we have κ u (δ) = u(ξ 0 ) -u(x 0 ) ≤ L 3 κ(δ).

2) If x 0 ∈ Ω and dist(x 0 , ∂Ω) ≤ Bδ, we can choose y 0 ∈ ∂Ω such that |y 0 -x 0 | = dist(x 0 , ∂Ω) ≤ Bδ. Then |ξ 0 -y 0 | ≤ (B + 1)δ ≤ ε 0 and |y 0 -x 0 | ≤ ε 0 . Since u is κ-continuous near the boundary , taking δ 3 small enough, it follows that |u(x 0 ) -u(y 0 )| ≤ C 0 κ(Bδ) and |u(ξ 0 ) -u(y 0 )| ≤ C 0 κ((B + 1)δ).

This implies that

κ u (δ) = u(ξ 0 ) -u(x 0 ) ≤ u(ξ 0 ) -u(y 0 ) + u(y 0 ) -u(x 0 ) ≤ C 0 κ(Bδ) + C 0 κ((B + 1)δ,
which by subadditivity implies κ u (δ) ≤ 2(B + 2)C 0 κ(δ). This proves the theorem.

Question : Is the main theorem true for for any modulus of continuity?

3.3. The case of quasi-plurisubharmonic functions. Let (X, ω) be a compact Hermitian manifold of complex dimension n. Let d be the geodesic distance on X associated to the metric ω. Let ϕ be an ωplurisubharmonic function on X. We define the modulus of continuity of ϕ as follows. For δ > 0 set

(3.11) κ ϕ (δ) := sup{ϕ(x) -ϕ(y) ; (x, y) ∈ X 2 , d(x, y) ≤ δ}.
On the other hand, we can define the local regularization R δ ϕ of ϕ on a neighbourhood of each point x 0 using a local chart (U, F ) centered at x 0 as follows : if F : U -→ C n is a biholomorphism from a neighbourhood U of x 0 to a bounded domain Ω

C n such that F (x 0 ) = 0. Then we define for x ∈ U

δ := F -1 (Ω δ ), R δ ϕ(x) := (ϕ • F -1 ) δ • F (x),
where Ω δ := {z ∈ Ω ; dist(z, ∂Ω) > δ} and (ϕ • F -1 ) δ is the standard regularization of the quasi-psh function ϕ • F -1 on Ω.

We consider a modulus of continuity κ : [0, l] -→ R + satisfying the following growth condition (3.12) ∃A > 0, lim sup

t→0 + 4nκ(At Aκ(t) < 1.
The following result was used in [START_REF] Demailly | Hölder Continuous Solutions to Monge-Ampère Equations[END_REF] for a Hölder modulus of continuity and was proved recently in [START_REF] Lu | Stability and Hölder regularity of solutions to complex Monge-Ampère equations on compact hermitian manifolds[END_REF] in that case. We will prove the following more general version using our previous results.

Theorem 3.4. Let κ be a modulus of continuity satisfying the condition (3.12) and ϕ : X -→ R be a bounded ω-plurisubharmonic function on X. Assume that for any point x 0 ∈ X there exists a local chart (U, z) contered at x 0 and constants C 1 > 0 and 0 < δ 1 << 1 such that for any x ∈ U and any 0 < δ < δ 1 ,

(3.13) R δ ϕ(x) -ϕ(x) ≤ C 1 κ(δ),
where R δ ϕ is a local regularization of ϕ in the local chart U . Then there exists a constant C 2 > 0 such that for any x ∈ X and y ∈ X, we have

|u(x) -u(y)| ≤ C 2 κ(d(x, y)).
Proof. It follows from our hypothesis that the function u is continous on X. Indeed fix an arbitrary point x 0 ∈ X and localize the problem in a neighbourhood of x 0 so that the inequailty (3.13) is staisfied. Fix a neighbourhood Y of x 0 in U and a biholomorphism F : Y -→ F (Y ) ⊂ C n from Y to a neighbourhood of the closed euclidean unit ball B ⊂ C n so that the point x 0 is sent onto 0. Since ω > 0 there exists a constant C > 0 such that β ≤ ω ≤ Cβ on Y , where F * (β)) is a multiple of the standard Kähler form on C n . Since u is ω-plurisubharmonic, and ω ≥ β, u is β-plurisubhramonic on Y and we can choose a local smooth potential w for β on Y so that the function v := ϕ + w is plurisubharmonic on Y . Since d ω (x, y) ∼ d β (x, y) on Y and R δ v = R δ ϕ + O(δ) on Y , we are then reduced to the case where v is a plurisubharmonic function on a neighbourhood of B which satisfies R δ v -v ≤ L 1 κ(δ) on a neighbourhood of B for some constant L 1 > 0. This implies that v is continuous on B, hence u is continuous on a neighbourrhood of x 0 . This proves the continuity of ϕ on X since x 0 was arbitrary.

In remains to prove the estimate on the modulus of continuity of ϕ. By continuity of u there exists (x δ , y δ ) ∈ X 2 such that d(x δ , y δ ) ≤ δ and κ ϕ (δ) = ϕ(x δ ) -ϕ(y δ ).

We want to show that there exists a constant C 2 > 0 such that lim sup

δ→0 + κ ϕ (δ) κ(δ) ≤ C 2 .
It's enough to show that the limsup along any sequence δ j → 0 is uniformly bounded by the same constant C 2 . Up to extracting a subsequence, we can assume by compactness that there exists a point x 0 ∈ X such that (x δ j , y δ j ) → (x 0 , x 0 ) as j → +∞.

Since R δ v -v = R δ ϕ -ϕ + O(δ), applying the previous localization process at the point x 0 , we are reduced to the case where v is a plurisubharmonic function on a neighbourhood Ω of B which satisfies R δ v -v ≤ C 1 κ(δ) on a neighbourhood of B. Then we can assume that all the points x δ j and y δ j belong to B for any j ∈ N.

From Lemma 3.3 there exists a constant L 2 > 0 such that for j > 1 large enough, we have v(x δ j ) -v(y δ j ) ≤ L 2 κ(δ j ), which implies that κ ϕ (δ j ) = ϕ(x δ j ) -ϕ(y δ j ) ≤ C 2 κ(δ j ). The theorem is proved.

Aknowledgements: We thank Chinh H. Lu for a careful reading of the preliminary version of this note and for useful comments on the proof of Lemma 3.3. We also thank Vincent Guedj and Henri Guenancia for useful discussions on the subject of this note.

  + δξ)dλ B (ξ), for x ∈ Ω δ := {x ∈ Ω ; dist(x, ∂Ω) > δ}, where λ B is the normalized Lebesgue measure on the euclidean unit ball B ⊂ R n .

S

  |u(x + δy) -u(x)| ≤ δ 0 S |Du(x + ty) • y|dσ(y)dt.

ANR project GRACK.