BOUNDARY-LAYER METEOROLOGY

(1)

BOUNDARY-LAYER METEOROLOGY

Han van Dop, September 2008

D 08-01

(2)

INTRODUCTION

This reader acts as a follow-up of the lecture notes of the Bachelor Hydro- dynamics and Turbulence course. Some of that material directly relating to boundary-layer meteorology, will be shortly discussed here.

1.1 Atmospheric thermodynamics

First Law of Thermodynamics

The change of energy dQ per unit mass of a closed system equals the sum of the change of the internal energy dU , and the amount of work, pdV :

dQ = dU + p dV. (1.1)

The specific heat at constant volume, or constant pressure, cv and cp are:

(dQ)v= (dU )v ≡ cv dT (dQ)p≡ cp dT

where (. . . )v denotes that volume is held constant, and (. . . )p indicates a constant pressure (cv,p in J kg⁻¹K⁻¹).

Equation of state

Again for a unit mass we have

pV = RT (1.2)

or, using ρ = 1/V :

p = ρRT. (1.3)

(4)

An alternative for the First Law can now be written as:

dQ = cp dT −1

ρdp. (1.4)

where cv+ R = cp, and R is the specific gas constant for dry air.

R = R^?/ma, with R^?, the universal gas constant (8.314 J K⁻¹ mole⁻¹) and mathe molecular weight of (dry) air (0.0288 kg). Thus R equals 288 J kg⁻¹K⁻¹.

Adiabatic lapse rate

The heat exchange between the atmosphere and its surroundings, for example by molecular conduction or by radiation, is often negligible, implying dQ = 0.

From the first law of thermodynamics and the hydrostatic equation, it follows that

cp dT + g dz = 0

−dT

dz = g

cp

≡ Γd≈ 10⁻²Km⁻¹ (1.5)

This vertical temperature gradient is commonly referred to as the (dry) adiabatic lapse rate.

Entropy

The entropy S is defined as

dS ≡ dQ T = cp

dT T − Rdp

p

= cp d ln T − R d ln p.

Upon integration, this yields:

S = cp ln T − R ln p. (1.6)

Potential temperature

In adiabatic processes, dS = 0. In a process (1 → 2), this means that :

(5)

cp ln T1− R ln p1= cp ln T2− R ln p2. (1.7) If we assign a constant pressure ps(e.g. 1000 mb) and temperature θ (potential temperature) to the reference level 2, then (1.7) can be written as:

c_p ln T − R ln p = c_p ln θ − R ln p_s, (1.8) where the subscript 1 is omitted. This equation can be rewritten to obtain

θ = T (ps

p)^κ, (1.9)

where κ = R/cp. This expression allows us to determine the potential temperature of any air mass as long as its temperature and pressure are known. In this way, a vertical profile of potential temperature can be constructed. Using (1.6), the entropy now becomes

S = c_p ln θ − R ln p_s, (1.10)

leading to the conclusion that, in case of an adiabatic process, there is a direct connection between the potential temperature, θ, and the entropy S. Further- more, Eq. (1.8) implies

dθ dz ≈ dT

dz + Γd, (1.11)

where the dry adiabatic lapse rate is given by

Γd= θ T

g cp

≈ g cp

. (1.12)

The above-mentioned equation poses a simple relation between the vertical temperature gradient dT /dz and the vertical gradient of the potential temperature dθ/dz.

Standard atmosphere

Using (1.5), an atmospheric temperature profile can be described:

T T0

= 1 − z

H, (1.13)

(6)

T0 being the surface temperature, and H = cpT0/g. We can now invoke the equation of state (p = ρRT ) and the hydrostatic equation to derive equations for vertical atmospheric profiles of pressure and air density:

p p0

= (1 − z

H)^cp^R (1.14)

ρ

ρ₀ = (1 − z

H)^cp^R⁻¹. (1.15)

The atmospheric scale height is denoted by H. This approach implicates that the atmosphere has a finite thickness (H ≈ 30 km) which is obviously not correct but nevertheless yields a reasonable estimate of the thickness of the atmosphere.

1.2 Statistical aspects of fluid mechanics

The length scales relevant in a turbulent flow cover a wide range. The macroscale and the smallest turbulent scale, known as the Kolmogorov microscale, are related through:

l/lK= Re³⁴ (1.16)

In the atmospheric boundary layer, typical values of l = 1000 m and lK= 0.001 m, yield a Reynolds number of O(10⁸). Solving all length scales up to 0.001 m on a domain of 10 km x 10 km x 1 km using a numerical program would require 10²⁰ grid points, whereas computation on a grid with 10¹⁰ points is currently feasible. Thus, turbulence equations cannot be solved numerically in the entire domain of length scales. As an alternative, turbulence equations could be solved for averaged quantities only. In this way, the range of length scales, and thus the required amount of grid points, can be reduced dramatically.

Stochastic variables, moments, probability density functions

In a turbulent flow, kinematic processes occur over a wide range of scales. As a first step in a statistical treatment of turbulence, the average of a variable (the first moment) can be determined. More detailed information of the flow can be obtained by calculating higher moments of variables, or combinations of variables (correlations). When analysing a turbulent flow in a statistical way, correlations between flow variables at different locations or at different times play a very important role. Examples of such correlations are the variance

u²≡ u(xxx, t)u(xxx, t) (1.17)

and the covariance

(7)

uw ≡ u(xxx, t)w(xxx, t), (1.18) being statistical properties of flow velocity at location xxx and time t. The overbar is used for the so-called ensemble-average of a quantity, i.e. the average obtained by taking the average of a large amount of realisations under the same boundary conditions.

If we treat velocity, pressure and other quantities in a turbulent flow as stochastic variables, it is possible to define the flow by its moments. Let u be any variable, and p(u) its probability density function. The moments are then given by

u^m≡ Z +∞

−∞

p(u)u^mdu (1.19)

A stationary stochastic process is defined as a process where the density function p is invariant under a time translation. If a stochastic process is stationary, p(u1, u2; t1, t2) = p(u1, u2; t1− t2) applies. The autocorrelation reads

u(xxx, t1) u(xxx, t2) = ρ(xxx, t1− t2).

Analogously, the spatial autocorrelation is defined as u(xxx1, t)u(xxx2, t).

In a homogeneous flow this function depends only on the separation:

ρ_s(xxx₁− xxx₂, t) = u(xxx₁, t)u(xxx₂, t).

where ρsis the (spatial) velocity autocorrelation function.

Reynolds decomposition and averaging processes

Reynolds decomposition

We can split the motion of a turbulent flow into a large scale and a small scale component (here denoted with a prime). In this way, the average of a variable and its fluctuation are discerned. Let uuu be the velocity. It can be split as follows

uuu = uuu + uuu⁰ (1.20)

’Averaging’, denoted by the overbar, means ensemble-averaging in this context.

We assume that the averaging process is a ’Reynolds operator’ which means that

(8)

u = u (and thus u⁰= 0), and that cu = cu (c is a constant) Averaging

Features of ensemble-averaging are:

• Ensemble-averaging removes all turbulence from the description of the flow: the variables are no longer chaotic. Any information about the structure of the turbulence vanishes.

• Transport by eddies has to be parameterised, i.e. only the statistical (average) properties of the turbulence are featured in the resulting equations.

Such a parameterisation requires a good knowledge of the flow.

• Averages vary much less in time and space than the turbulence itself. Grid point distances can therefore be chosen relatively large.

Time-averaging

The average value of a velocity variable over a period T is defined as

UUUT(xxx, t) = 1 T

Z t+¹₂T t−¹₂T

uu

u(xxx, t⁰)dt⁰, (1.21)

If T incorporates all turbulence time scales, the time-average would be a good approximation of the ensemble-average (U will then be independent of T ).

Volume-averaging

Analogous to time-averaging, a variable can also be averaged over a certain volume. The local values are split into a spatial average and a deviation from this average

a = {a} + a⁰ (1.22)

where the volume-average is defined as

{a} ≡ 1

∆V Z

∆V

a dV. (1.23)

In more general terms, the volume-average can also be expressed as

{a}(x) ≡ Z +∞

−∞

G(x⁰− x)a(x⁰)dx⁰ (1.24)

(9)

where the filter function G is given by, for example,

G(x − x⁰) = 1 if |x⁰− x| < L

= 0 otherwise.

Important properties of (1.23), that differ from the properties of ensemble- averages, are

{{a}{b}} 6= {a}{b}

{{a}} 6= {a}

(1.25) In practical applications these inequalities are often ignored, since errors made by assuming that these quantities are equal appear to be small. Characteristics of volume-averaging are:

• Volume-averaging only removes the scales smaller than V¹³; it conserves the larger scales.

• A so-called closure is needed for the scales smaller than V¹³. The advantage of volume-averaging, compared to ensemble-averaging, is that the closure hypothesis is less critical, since the small scales do not carry much energy and have relatively well-known properties.

(10)

Chapter 2

CONSERVATION LAWS

We summarize the basic equations as follows:

The continuity equation

dρ dt + ρ∂ui

∂xi

= 0. (2.1)

The momentum equation

ρ∂ui

∂t + ρuj

∂ui

∂xj

= −∂p

∂xi

− ρgδi3+ µ ∂²ui

∂xj∂xj

. (2.2)

The temperature equation

∂θ

∂t + uj

∂

∂xj

θ ≡ dθ

dt = κ ∂²θ

∂xj∂xj

. (2.3)

2.1 Governing equations

As a good approximation, the continuity equation can be expressed as

∂ui

∂x_i = 0, also implying that

(11)

1 ρ

dρ dt ≈ 0.

The density along the fluid-particle trajectories is approximately constant. The flow can therefore be considered incompressible.

In an adiabatic process, dρ

dt = 1 c²

dp dt

applies, where c (the velocity of sound) stands for (^c_c^p

v p

ρ)¹². When combined with the Bernoulli equations, this yields

∂ui

∂xi

= 1 c²

(d¹₂ ui2

dt + gu3

) .

Choosing U and L as scale sizes of uiand xirespectively, then

L U

∂ui

∂xi

= O(u²

c²) + O(g L c² ).

Using the Boussinesq approximations, the right-hand side can be neglected.

Boussinesq approximations

The above approximations, u/c << 1 and L << c²/g, which imply that atmospheric flow can be considered incompressible, form part of the so-called Boussinesq approximations. This however does not mean that density differences are dynamically unimportant. They certainly are in combination with gravity. The resulting equations are referred to as the Boussinesq equations.

They originate from a number of subtle arguments which will be given below.

We will suppose a reference state of the atmosphere (denoted by the subscript 0). The atmosphere is at rest (u_i= 0), and horizontally homogeneous regarding temperature, density and pressure. We further assume that thermodynamical deviations from the reference state in a realistic atmosphere are small. We further suppose that the actual atmospheric state does not deviate much from the reference state, viz.

p −→ p0+ ˜p ρ −→ ρ₀+ ˜ρ θ −→ θ0+ ˜θ ui−→ 0 + ˜ui,

where ˜p, ˜ρ, ˜θ and ˜ui denote small deviations from the reference state. When we substitute this in the equation of state, p = ρRθ, (we assume that in the ABL

(12)

θ ≈ T , see Eq. (1.9)). We obtain in zero order:

p0= ρ0Rθ0, (2.4)

and to first order

˜

p = R(ρ₀θ + ˜˜ ρθ₀).

We assume, confirmed by observations and given the approximations already made, that

˜

p << ˜θ, ˜ρ, so that

˜ ρ ρ0

∼ − θ˜ θ0

. (2.5)

Now we make the same substitutions in (2.2):

(ρ0+ ˜ρ)d ˜ui

dt = −∂p0

∂xi

− ∂ ˜p

∂xi

− (ρ0+ ˜ρ)gδi3+ µ∂²u˜i

∂x²_j . In zero order this yields

∂p0

∂xi

= −ρ0 g δi3, (2.6)

the hydrostatic equation, and to first order:

ρ0

d ˜ui

dt = − ∂ ˜p

∂xi

− ˜ρg δi3+ µ∂²u˜i

∂x²_j , which can be rewritten, using (2.5), as

d ˜ui

dt = − 1 ρ₀

∂ ˜p

∂x_i + θ˜

θ₀g δi3+ ν∂²u˜i

∂x²_j ,

where ν = µ/ρ₀, the kinematic viscosity. Back substitution of the original variable θ yields

d ˜ui

dt = − 1 ρ0

∂ ˜p

∂xi

+ (θ − θ0

θ0

)gδi3+ ν∂²u˜i

∂x²_j . (2.7)

This is a familiar expression of the Boussineq equations for a shallow boundary layer. Essential is that density (temperature) deviations are only important in combination with gravity. Otherwise density can be considered constant. So the Boussinesq-approximations for the (shallow) ABL include:

(13)

u c << 1 L << c²

g θ ≈ T

˜ p p<< 1

and we summarize the equations in Boussinesq form:

∂ ˜ui

∂x_i = 0. (2.8)

The momentum equation

∂ ˜ui

∂t + uj

∂ ˜ui

∂x_j = − 1 ρ₀

∂ ˜p

∂x_i +(θ − θ0)

θ₀ gδi3− 2ijkΩjuk+ ν ∂²u˜i

∂x_j∂x_j. (2.9) For the sake of completeness, the Coriolis term has been added, making the equations apply for a rotating coordinate system, with an angular velocity ΩΩΩ (in rad/s).

The temperature equation

The temperature equation, neglecting molecular conduction, is changed to d˜θ

dt = 0. (2.10)

Together, eqs.(2.8, 2.9 and 2.10) constitute a set of 5 equations (with unknown variables ˜u_i, θ and ˜p) that serves as a starting point for the study of the atmospheric boundary layer.

(14)

Chapter 3

THE REYNOLDS EQUATIONS

We will suppose that a flow consists of a laminar, average flow, and, superim- posed, turbulent fluctuations, the so-called Reynolds decomposition (see section 1.2).

3.1 The average flow

We will adopt the Navier-Stokes equation, Eq. 2.9,

∂ ˜ui

∂t + ˜uj

∂ ˜ui

∂xj

= − 1 ρ0

∂ ˜p

∂xi

+(θ − θ0) θ0

gδi3− 2ijkΩjuk+ ν ∂²u˜i

∂xj∂xj

. (3.1) and the continuity equation

∂ ˜u_i

∂xi

= 0. (3.2)

We decompose the velocity ˜u_i, the pressure ˜p and the temperature θ in a mean and fluctuating (turbulent) component. After substitution of ˜ui = Ui+ ui,

˜

p = Π + π and θ = Θ + θ (with this substitution we redefine from hereon θ as a temperature fluctuation) and averaging, we obtain

∂U_i

∂t +U_j∂U_i

∂xj

= −1 ρ0

∂Π

∂xi

+(Θ − θ₀) θ0

gδ_i3−2_ijkΩ_jU_k+ν ∂²U_i

∂xj∂xj

−∂u_iu_j

∂xj

, (3.3) and

∂Ui

∂xi

= 0. (3.4)

(15)

Figure 3.1: Frictional flow along a flat surface

The extra terms, uiuj, which represent turbulent momentum fluxes, are called Reynolds stresses and are expressed in known variables. This procedure is called closure (or ’K-theory’) and will be extensively treated in section 3.5:

u_iu_j = −K_m ∂U_i

∂xj

+∂U_j

∂xi

(3.5)

The proportionality constant K_mtakes the dimension of viscosity (m²s⁻¹), but is entirely determined by the turbulence itself.

Let l and u be the length and velocity scales of the turbulence, then it seems reasonable to suppose that Km∝ lu. This can be nicely illustrated by examining a flow along a flat surface.

Neutral boundary-layer flow

The average flow is stationary and homogeneous in the x-direction (^∂U_∂t,^∂U_∂x = 0).

The equations of motion (meglecting Coriolis forces) are (see 3.3)

∂

∂z ν∂U

∂z − uw

= 1

ρ0

∂Π

∂x (3.6)

1 ρ₀

∂Π

∂z −∂w²

∂z + g = 0. (3.7)

By differentiating (3.7) to x, we find that _ρ¹

0

∂Π

∂x cannot be a function of z. Now, integrating the first equation yields:

ν∂U

∂z − uw

∝ z. (3.8)

(16)

Suppose that h is the ”thickness” of the boundary layer (the height at which

∂U/∂z ≈ 0)). Consider two cases:

a. Laminar flow close to the wall (uw = 0).

Eq. (3.8) gives ν∂U

∂z = az + b, (3.9)

where the constants a and b follow from the boundary conditions at z = 0 and z = h. At z = h we assume ∂U/∂z = 0. b represents the friction force at the surface (with dimension velocity squared). This defines the ’friction velocity’ u_? as

lim

z→0ν∂U

∂z ≡ u²_?, (3.10)

so that we finally get ν∂U

∂z ∝ u²_?(1 − z

h), (3.11)

implying that the velocity profile is linear for z << h. In this (thin) layer molecular friction dominates and thus Re = O(1). If the thickness of that layer is δ we have thus

u_?δ ν ≈ 1,

which yields δ ≈ ν/u_? for the thickness of the ’laminar sub-layer’.

b. Turbulent flow: the logarithmic wind profile

Above the laminar sub-layer (ν/u_?. z . h) turbulent friction dominates:

|uw| >> ν^∂U_∂z. Eq. 3.8 now reads

−uw ∝ z

which, applying the boundary conditions −uw = 0 at z = h and −uw = u²_? at z = ν/u?, yields

−uw ∼ u²_?(1 −z

h). (3.12)

(note that ν/u_? << h). If z/h << 1, the vertical transport of momentum (−uw ≈ u²_?) is approximately constant. This defines another sub-layer, the

’constant flux layer’, ν/u_?. z << z/h.

(17)

We can now use the closure relation (3.5) to express the Reynolds stress in terms of the average gradient in this layer:

−uw = Km

∂U

∂z = u²_?. (3.13)

According to (3.13), the velocity profile would be linear (using a constant value for Km). This has not been confirmed experimentally: instead, logarithmic profiles are found. This only follows from (3.13) if the ’eddy’ diffusion coefficient K is proportional to z. If we define

Km= κzu?,

where κ is the Von Karman constant), the solution of (3.13) will be U

u?

= 1 κln z

z0

, (3.14)

This is a logarithmic profile¹, where z0 is an integration constant. The proportionality constant κ has a value of ≈ 0.4. The quantity z0 depends on the roughness of the wall, and is called the roughness length. For an aerodynamically smooth wall, z0 is defined as the height at which the corresponding Reynolds number has the value of 1:

(Re)z0≡ u_?z₀

ν , (3.15)

leading to z0 = ν/u?. If an aerodynamically rough wall is characterised by irregularities of height h, the Reynolds number is given by

(Re)h≡ u?h

ν . (3.16)

A surface is called smooth when Reh < 1 and rough when Reh > 1. A couple of typical values for z0 are listed in table 3.1. The roughness of a rough water surface is treated later in these notes.

Under neutral conditions and over a large area of varying surface types, the wind profile is approximately logarithmic from a height z0to a couple of hundreds of metres.

If the area is covered with higher objects (trees), it is possible to define a new surface where the wind speed is 0. The wind profile is then given by

U u?

= 1

κlnz − d z0

, (3.17)

1There is as yet no exact derivation for the logarithmic wind profile. In section 4.3, it is made plausible that logarithmic profiles occur in boundary layers

(18)

Table 3.1: Typical values for the roughness length z₀ Surface type Roughness length (m) Smooth water/ice 10⁻⁴

Short grass 10⁻² Low vegetation 0.05 Countryside 0.20 Low built-up area 0.6 Forests/cities 1 − 5

where d is the so-called displacement height. The displacement height equals roughly 80% of the object height (see figure 3.2). Equation (3.17) is very suitable to calculate the wind speed at a given height, if the wind speed is known on any other height:

U₂ U1

=ln(z₂− d)/z0

ln(z1− d)/z0

.

Drag coefficient

If we square (3.14) and subsequently multiply it by the density, we will find the surface shear stress:

τ ≡ ρ u²_?= ρ κ²

ln²(z/z0) U². The drag coefficient is thus defined as

Cd= κ²

ln²(z/z0), (3.18)

implying that

τ = ρ C_d U², (3.19)

This is a simple relationship to estimate the surface shear stress from the average wind speed under neutral circumstances.

Power law

In practice, an algebraic formula for the wind profile is often applied:

(19)

Figure 3.2: Sketch of the wind speed profile over a homogeneous forest.

U2

U₁ = z2

z₁

^p

. (3.20)

In neutral conditions, the following rule approximately holds:

p ≈ ln⁻¹ √z1z2

z₀

. (3.21)

This relation can be derived by the requirement that at the geometric mean value of z₁and z₂,p(z₁z₂), the wind velocity and the first derivative are equal in both formulations.

A frequently used, but not necessarily correct value of p is 1/7 (based on z1 ∼ z2∼ 10m and a roughness length of 8 cm).

3.2 The mean-flow energy equation

The starting point is (3.3). We shall neglect the Coriolis force since it does not play a role in energy budget considerations:

∂Ui

∂t + Uj

∂Ui

∂xj

= − 1 ρ0

∂Π

∂xi

+(Θ − θ0) θ0

gδi3+ ν ∂²Ui

∂xj∂xj

−∂uiuj

∂xj

, Multiplying by U_i yields:

(20)

∂E

∂t + Uj

∂E

∂x_j

= − 1 ρ0

Ui

∂Π

∂xi

+Ui(Θ − θ0) θ0

gδi3+ ∂

∂xj

{νUi

∂Ui

∂xj

− Uiuiuj} + uiuj

∂Ui

∂xj

− ν(∂Ui

∂xj

)², (3.22) where E ≡ ¹₂U_i². If the Reynolds number is high, the viscous terms can be omitted.

When integrated over the volume of the flow, the result is dEtot

dt = Z

V

uiuj

∂Ui

∂xj

dV + Z

V

Ui(Θ − θ0) θ0

gδi3 dV. (3.23)

The integrand of the first integral is negative in a shear flow, and is called the deformation work. This term provides the translation of the average flow energy to its fluctuations (i.e. the coupling between average flow and turbulence). The second integral represents the conversion of potential into kinetic energy by the mean motion in a stratified flow.

3.3 The turbulent kinetic energy equation

The energy transfer from the average flow to turbulence is provided by the deformation work. We will now derive an expression for the average kinetic energy of the fluctuations, q ≡ ¹₂ui2. We shall start from (2.9):

∂ ˜u_i

∂t + ˜u_j∂ ˜u_i

∂xj

= − 1 ρ0

∂ ˜p

∂xi

+(θ − θ₀) θ0

gδ_i3+ ν ∂²u˜_i

∂xj∂xj

, and for the temperature

∂θ

∂t + ˜u_j ∂θ

∂x_j = κ ∂²θ

∂x_j∂x_j,

As in the derivation of the equation for the mean energy (3.3), we make again the substitutions ˜ui = Ui+ ui, etc. In order to arrive at an equation for the fluctuations, we substract from this equation the equation for the mean flow (3.3). We multiply the resulting equation with ui. Taking the average of this equation, we finally obtain the equation for the energy of the fluctuation q ≡

1 2ui2:

∂q

∂t+Uj

∂q

∂xj

= −uiuj

∂Ui

∂xj

+ g θ0

wθ − ∂

∂xj

{¹₂uiuiuj+1 ρ0

πuj−ν∂¹₂u²_i

∂xj

}−. (3.24)

(21)

The first term on the right-hand side of (3.24) is the mechanical production term, supplied by the average flow (normally positive). The second term represents the thermal production of turbulent fluctuations due to differences in density of the flow. In a gravity field, this term can transform potential energy into kinetic and vice versa. The ratio between mechanical and thermal production of turbulence is defined as the flux-Richardson number:

Rif =

g θ₀wθ u_iu_j^∂U_∂xⁱ

j

(3.25)

We discern three cases:

wθ > 0 unstable flow Ri_f < 0 wθ = 0 neutral flow Rif = 0 wθ < 0 stable flow Rif > 0

The use of the flux-Richardson number can be illustrated by the atmospheric boundary layer that is heated by the Earth surface during the day, so that wθ > 0. At daytime, the boundary layer is thus generally unstable. During the night, the Earth surface will cool down by radiating out its energy. Thus, wθ < 0, meaning that the boundary layer will be stable with little turbulence.

We also distinguish the gradient Richardson number, Ri, by applying K-theory to the fluxes in Eq. 3.25. Similar to Eq. 3.5 we may write for wθ:

wθ = −Kh

∂Θ

∂z, (3.26)

so that Eq. 3.25 becomes Ri_f = K_h

Km

Ri, (3.27)

where Ri is given by

Ri =

g θ0∂Θ/∂z (∂Ui/∂xj+ ∂Uj/∂xi)^∂U_∂xⁱ

j

. (3.28)

Apart from the flux-Richardson number, we also use the bulk-Richardson number, which follows from (3.28), by approximating gradients by finite differences and assuming a uniform windfield, U (z),

Rib= g θ₀

∆Θ

(∆U )²∆z. (3.29)

(22)

Returning to Eq. 3.24, we observe that the redistribution term between braces contains 3 energy fluxes: by turbulent fluctuations, by pressure fluctuations, and by molecular fluctuations (the latter has already been omitted in the equation).

The last term represents the energy dissipation. It is the only negative term, and therefore necessarily of the same order of magnitude as the other terms. If not, then q could grow unlimitedly.

In a neutral, semi-stationary situation, (3.24) implies

P ≡ −uiuj

∂Ui

∂xj

≈ ν(∂ui

∂xj

)²≡ ,

so P ≈ . The production and dissipation of turbulent energy are approximately balanced. Also note that P = O(^u_l³), which means that

= O(u³

l ). (3.30)

3.4 Spectra

Turbulence can be viewed upon as a set of eddies of different sizes between l_K and L. In a turbulent flow, energy of a particular scale is transferred towards larger and smaller scales. The average flow provides the energy for the turbulence at the macroscale. The larger eddies are unstable and disintegrate into smaller eddies of various sizes. This process is repeated several times. The smallest eddies will ultimately lose its energy by molecular viscosity. The net effect is that eddies tranfer energy from larger to smaller scales, known as the energy cascade.

A turbulent velocity field consists of a superposition of fluctuations with a wide range of temporal and spatial scales. To analyse the contribution of each scale, a Fourier analysis can be exploited. Let the autocorrelation of a continuous velocity function u(t) be defined as

ρ(τ ) = u(t) u(t + τ ) = lim

T →∞

1 T

Z +^T₂

−^T₂

u(t) u(t + τ ) dt. (3.31) We have assumed stationary conditions so that ρ depends on the time difference τ only. The energy spectrum S(ω) is, by definition, the Fourier transform of the correlation function:

S(ω) ≡ 1 2π

Z +∞

−∞

e^−iωτ ρ(τ ) dτ, (3.32)

with the associated inverse operation ρ(τ ) ≡

Z +∞

−∞

e^iωτS(ω) dω, (3.33)

(23)

From (3.33), it immediately follows that for τ = 0 u²≡R+∞

−∞ S(ω) dω.

The lefthand side represents the turbulent kinetic energy and this equation shows why S is actually called the energy spectrum, since S(ω) equals the contribution to u²of S between the frequencies ω and ω + dω.

(24)

In practice the energy or power spectrum of the function u(t) is determined in a more direct way: we rewrite Eq. 3.31 as

ρ(τ ) = lim

T →∞

1 T

Z +∞

−∞

v(t) v(t + τ ) dt, (3.34)

where v(t) is defined as

v(t) = u(t) for|t| < T /2 v(t) = 0 otherwise.

The Fourier tranform of v(t) is

˜ v(ω) = 1

2π Z+∞

−∞

e^−iωt v(t) dt, (3.35)

with the inverse transformation v(t) =

Z +∞

−∞

e^iωtv(ω) dω.˜

We use the last expression to rewrite Eq. 3.34 as ρ(τ ) = lim

T →∞

1 T

Z +∞

−∞

v(t) Z+∞

−∞

˜

v(ω) e^{iω(t+τ )}dω dt.

Rearranging terms we get ρ(τ ) = lim

T →∞

1 T

Z +∞

−∞

˜ v(ω) e^iωτ

Z+∞

−∞

v(t)e^iωtdt dω, or

ρ(τ ) = lim

T →∞

2π T

Z +∞

−∞

˜

v(ω) ˜v(−ω) e^iωτ dω.

Since ˜v(−ω) = ˜v^?(ω) we have

ρ(τ ) = lim

T →∞

2π T

Z +∞

−∞

|˜v(ω)|²e^iωτ dω.

Comparing this result with Eq. 3.33 yields S(ω) = lim

T →∞

2π T |˜v(ω)|², which can be rewritten as

S(ω) = lim

T →∞

2π T

˛

˛ Z+^T₂

−^T 2

u(t) e^iωtdt

˛

2

. (3.36)

This equation, known as the Wiener-Khintchin theorem, provides a direct relation between the velocity signal u(t) and its power spectrum. The Fast Fourier Transform (FFT) is an efficient numerical way to determine the spectrum of an arbitrary stationary process based on a discrete representation of Eq. 3.36.

Some of the properties of S are

• S(ω) = S(−ω)

(25)

• S(ω) is real and ≥ 0.

• It also follows from the definition that S(0) = _π¹R∞

0 ρ(τ ) dτ = _π¹(u² Tu), where Tu denotes the time scale of the signal u(t).

Spectral gap

During the 1950s, the ideas about boundary-layer meteorology were dominated by the assumption that (small-scale) turbulent kinetic energy mainly occurred at scales comparable to the boundary layer height (1-2 km). At larger scales, hardly any energy would be available, whereas the energy would increase at mesoscale and sub-synoptic scales. This implied that the turbulence energy spectrum would contain a minimum, separating boundary-layer turbulence from turbulence at larger scales. This would also justify the process of Reynolds decomposition. Now that much more experimental data are available, this view has become challenged. An example of a spectral analysis of the energy, recorded during an airplane flight, is given in figure 3.4. In this figure, the spectral energy of the u− and v velocity components continue to increase at larger scales. Other turbulent variables, like temperature, water vapour content, and liquid water concentration show comparable behaviour. Only the spectrum of the vertical velocity w seems to level off towards larger scales.

That the spectral energy at larger scales does not vanish, has some inconvenient implications for the integral properties of the spectrum. For example, the total variance of u, being defined as

u²≡ Z ∞

0

Eu(k)dk, (3.37)

cannot be determined if Eu(k) continues to grow for k → 0. In these cases, the ogive is defined as

Og_u(k₀) ≡ Z ∞

k₀

E_u(k)dk, (3.38)

so that integral properties of the spectrum can be determined anyway. The choice of k₀ is an arbitrary one.

3.5 Closure

Expressing the stress in terms of the average flow is called closure. More generally, closure means that higher moments are expressed in terms of lower ones.

It is called closure because the set of equations of motion can be solved after closure. In the previous sections, some simple examples of 1^st-order closure have already been demonstrated. More elaborate ways of closure exist, which will be shortly discussed below.

(26)

Figure 3.3: Kaimal’s schematic representation of the energy spectrum. (A) production subrange; (B) the inertial subrange, where production and dissipation are not important; (C) dissipation subrange. The second Kolmogorov hypothesis conjectures that in the inertial subrange, the dissipation and the wavenumber k are the dominating parameters. The spectrum E(k) can be determined using dimensional analysis based on and k only. The result is E(k) = CK²³k⁻⁵³.

First-order closure (K-theory)

First-order closure is often applied since it is a fast method. It implies that the

0K-theory’ is used to determine the Reynolds-stress term (see, e.g., Eq. 3.5 ).

The diffusion coefficient K is typical for the flow and may thus change in time and space. At the surface layer, we acceptably showed that Km = κzu? (see chapter 3.3, Eq. 3.13). Above the surface layer, there are many possible ways to follow. At the top of the boundary layer, the vertical exchange is small, leaving Kmapproximatly zero. Somewhere within the boundary layer, Kmshould have a maximum. A possible formulation for Km is:

Km= l²

∂U

∂z

2

+ ∂V

∂z

2¹₂

. (3.39)

In this equation, l is the mixing length, that will be proportional to the height in the surface layer, l = κz, and converge to a constant value, say λ, for greater heights. This choice for K_mwill give the desired behaviour in the surface layer, where V v 0 and −uw are constant (u²?). It now follows from

uw = −Km

∂U

∂z (3.40)

(27)

Figure 3.4: 1-D spectra of the horizontal (u and v) and vertical (w) wind-speed components as a function of the wave number k. Measurements from an airplane at approx. 150 m above the sea.

(28)

after substitution in (3.39), that

∂U

∂z = u?

κz,

(see Eq. 3.13) which provides the desired logarithmic wind profile. The height- dependent growth of the length scale can be delimited by choosing the following parameterisation:

l = κz

1 + ^κz_λ (3.41)

where λ is an empirical heigth scale (typically ∼100 m). The eddie diffusion parameter Km (3.39) can easily be corrected for the stability by adding an empirical function F :

K = l²

∂U

∂z

2

+ ∂V

∂z

2¹₂

F (Ri), (3.42)

F being a function of the stability through the Richardson number. This parameterisation from 1979 is still used in present-day global circulation models like the ECMWF-model, that calculate weather forecasts.

1.5-order closure

This type of closure makes use of the energy equation (3.24), which is a second- moment equation. Suppose that Kmobeys

K_m≈ lq¹², (3.43)

where q is given by the turbulent kinetic-energy equation (3.24). The terms present in that equation can be approximated in the following way, if the flow is assumed to be horizontally homogeneous:

−uiuj

∂Ui

∂x_j = −Km

∂U

∂z

² g

θ₀wθ = −g θ₀ Kh

∂Θ

∂z

∂

∂xj

1

2uiuiuj+πuj

ρ0

= −∂

∂zKm

∂q

∂z

= q³²l⁻¹.

Substituting the above equations into the stationary-energy equation yields:

(29)

0 = K_m ∂U

∂z

² + g

θ0

K_h∂Θ

∂z + ∂

∂zK_m∂q

∂z − q³²l⁻¹. (3.44) An equation for the turbulent kinetic energy q emerges. The set of equations (3.3, 3.5, 3.41, 3.43 and 3.44) can now be solved. The name ’1.5-order closure’

is used, since only one second-order equation is used for the closure algorithm.

Second-order closure

In this case, all second moments, like u_iu_j, u_iθ, are incorporated in the closure algorithm. These equations contain terms of the third order, which are either parameterised or neglected. As a result, the number of equations that are to be solved increases rapidly. Sometimes, this method yields better results. On the other hand, the mathematical complexity, the large number of empirical constants and the lack of a solid theory make second-order closure techniques less favourable.

Large Eddy Simulation (LES)

LES makes use of a fundamentally different approach. The variables in the Navier-Stokes equation are averaged over a small volume. Contrary to the variables in the Reynolds-averaged equations, the LES-variables describe the turbulent behaviour of the flow at scales larger than the scales of the averaging process. The starting point is the momentum equation (see 3.1), where we take for simplicity a neutral flow and neglect the Coriolis-term and molecular viscosity. The volume-average is:

∂{ ˜ui}

∂t + { ˜uj}∂{ ˜ui}

∂xj

= − 1 ρ0

∂{˜p}

∂xi

+∂Rij

∂xj

, (3.45)

where Rij is a tensor that describes the impact of the flow behaviour within the volume on that outside the volume. This is called the subgrid scale stress, defined as Rij = { ˜ui}{ ˜uj} − { ˜uiu˜j}. Comparable to Reynolds-averaging, this term is approximated by relating it to volume-averaged variables:

Rij = −νt

∂{ ˜ui}

∂x_j +∂{ ˜uj}

∂x_i

(3.46) In a mathematical sense, LES is almost identical to the method of Reynolds- averaging. The difference between the two processes can be found in the way the eddy viscosity is formulated. Reynolds-averaging appoints an order of magnitude u.l to K, so that the Reynolds number, based on K, is of the order of 1.

The choice for the eddy viscosity in LES, νt, depends on the volume over which the averaging process takes place. This volume is usually of the same order of

(30)

magnitude as the grid point distance of the numerical scheme, ∆. The Reynolds number, ul/νt≈ l/∆ is now much larger than 1 (approximatly 10²to 10³using current computer capacity).

A well-known choice for the eddy viscosity is formulated by Smagorinsky (1963), which is based on a mixing hypothesis (like the K-theory):

ν_t= (C_s∆)²|S| (3.47)

where Csis the Smagorinsky constant (Cs= 0.2), |S| = (2SijSij)¹² and Sij= ¹₂_{∂{ ˜}_u

i}

∂xj +^{∂{ ˜}_∂x^u^j^}

i

, the deformation tensor.

The most important applications of LES can be found in the simulation of convective turbulence, where transports are dominated by large-scale motions in the atmosphere.

Direct Numerical Simulation (DNS)

This technique aims to find solutions for the Navier-Stokes equation itself, with- out introducing closure or averaging. At large Reynolds numbers, this approach faces difficulties, since the range of turbulence scales increases rapidly. In order to find numerical solutions, the grid points distances have to be approximately equal to the smallest scale. For calculating practical situations, a large number of grid points are required. Flows with Reynolds numbers up to 100-1000 can be simulated with currently available computer power.

(31)

Chapter 4

THE ATMOSPHERIC BOUNDARY LAYER

4.1 Phenomenology

Boundary-layer meteorology comprises the dynamics and physics of the atmospheric layer that is closest to the Earth surface. Exchange processes between the Earth surface and the ’free atmosphere’ occur in the boundary layer. The state of the boundary layer is influenced by flow in the free atmosphere on the one hand, and by boundary conditions imposed by the Earth surface on the other hand.

Dominating boundary layer processes are the vertical exchange of momentum τ = −ρ₀ uw, heat H = ρ₀ c_pwθ and water vapour E = ρ₀ wq.

The free atmosphere (see figure 4.3) is the layer above the uppermost boundary of the (turbulent) boundary layer.

The behaviour of the boundary layer has an enormous impact on values of the maximum and minimum temperature, wind speed gradient (wind shear), fog and cloudiness. In aviation and architecture (wind engineering), both wind shear, fog and the occurrence of wind gusts are important factors to keep in mind. Turbulence, radiation, surface properties and thermodynamics are the key ingredients of boundary-layer meteorology. We will first give a qualitative overview.

Slightly above the Earth surface, a thin laminar layer exists: the ’viscous sub- layer’. This layer adjusts the balance between the Earth surface and the surface layer. The surface layer (or inner layer) is a couple of tens of metres thick. The surface layer is topped by the (planetary) boundary layer, sometimes also referred to as the Ekman layer or outer layer. The transition between the surface layer and the Ekman layer is smooth. The scales in the surface layer are very different from those in the Ekman layer.

The boundary-layer dynamics are not only influenced by the average horizontal

(32)

Figure 4.1: Schematic overview of the troposphere.

(33)

flow, but also by turbulent processes. The most important of these processes are:

a. Mechanical turbulence

b. Buoyant or convective turbulence.

The stability of the atmosphere influences the boundary layer, which has a typical height of 100-2000 m. The flux-Richardson number is a measure for the atmospheric stability. It is the ratio between buoyancy production and mechanical production

Rif =

g θ0wθ uiuj∂Ui

∂x_j

(see section 3.3).

The surface layer is heated when the surface heat flux wθ0 is positive. This happens during daytime as a consequence of solar radiation. Convection may also take place above a warm sea surface, over which cooler air is flowing.

If the heat flux is negative, the boundary layer will cool down. This happens during the night. These two situations differ substantially, having quite some implications for the dynamics (e.g. the growth of the boundary layer height:

if wθ > 0, then Ri_f < 0, so the atmosphere is unstable. This implies that gravity is generating turbulence, and the boundary layer height increases as a consequence. On the other hand, if wθ < 0, so Ri_f > 0, the boundary layer is stable. Gravity will suppress turbulence. The boundary layer height will not change notably. We conclude that there is a daily cycle of boundary layer height over land surfaces (see figure 4.2).

• Unstable boundary layer

Incoming shortwave radiationwill heat the Earth surface which will heat and lift the surface-layer air: convection (thermals) will start to emerge.

In this manner, the boundary layer becomes unstable. Normally, there is a stable layer on top of the boundary layer, which limits the growth of convection (see figure 4.3). When convective cells penetrate into a stable layer, they will lose their kinetic energy and won’t ascend any further.

These convective cells mix the air of the boundary layer with the air above. In this way, the thickness of the boundary layer increases. This process is called entrainment. The strong turbulence makes the mixing effective, resulting in a rather uniform distribution of momentum, heat and water vapour.

A (moist) convective cell that reaches its LCL (Lifting Condensation Level) will condensate and form a cumulus cloud. As soon as the cell reaches its LFC (Level of Free Convection) the clouds can grow upward until the next stable atmospheric layer. The limit of this entrainment process is at the level of the tropopause, an extremely stable layer at a height of about 10-15 kms.

(34)

Figure 4.2: Typical diurnal progress of the boundary layer over a land surface.

(From: Stull, An Introduction to Boundary-Layer Meteorology, 1988).

Figure 4.3: Characteristic profiles of average potential temperature, wind speed, vapour and a random trace gas in an unstable boundary layer (From: Stull, An Introduction to Boundary-Layer Meteorology, 1988).

(35)

Figure 4.4: Typical profiles of potential temperature and wind speed in a stable boundary layer (From: Stull, An Introduction to Boundary-Layer Meteorology, 1988).

• Stable boundary layer

A stable boundary layer emerges when the net longwave radiation of the Earth surface is negative, which cools down the air above it. Vertical motions are often hampered in a stable boundary layer, since the (negative) buoyancy slows them down. The turbulent kinetic energy and the length scales are much smaller than in an unstable boundary layer. The turbulence is not well-structured. The height of a stable boundary layer is determined by the balance between turbulence production and dissipation.

Since turbulence is suppressed, little mixing takes place in the boundary layer. Gradients of heat, momentum, and water vapour are therefore much larger (see figure 4.4).

4.2 The surface layer; Monin-Obukhov theory

The Monin-Obukhov length

In a non-neutral boundary layer, the surface heat flux (wθ)₀ is not zero. To- gether with the variables z and u_?, the surface heat flux, or rather the buoyancy

g

θ0(wθ)₀plays a significant role. Since u_? and _θ^g

0(wθ)₀ can be more or less constant over an hour or longer, they are considered to have a big influence on the dynamical structure of the surface layer. Based on these variables, we can derive a stability parameter, the Monin-Obukhov (MO) length scale L:

L = − u?3

κ(g/θ₀)(wθ)₀.

The turbulent kinetic energy equation (3.24) can be simplified to

(36)

0 = −uw∂U

∂z + g θ0

wθ − ,

and it can be made dimensionless through multiplication by κ z/u³_?:

0 = −uw u²_?

κ z u?

∂U

∂z − z L − φ.

We used the definition of L to arrive at the latter equation, and denoted the dimensionless dissipation of energy by φ. Since in the surface layer −uw ≈ u²_?, we find:

0 = φm− z L− φ,

in which we wrote the dimensionless wind speed profile as

φ_m= κ z u?

∂U

∂z. (4.1)

In this simplified form, the energy equation can be interpreted as a balance between mechanical and buoyancy production on the one hand, and the dissipation of energy on the other. The equation is expressed in only a few scales (z, L, and u_?). From dimensional analysis, we can also define a dimensionless temperature gradient φ_h:

φh(z/L) =κz θ_?

dΘ

dz. (4.2)

The temperature scale θ_? is defined as θ_?≡ −wθ₀/u_?.

The functions φ_m and φ_h have been determined experimentally. In case of instability (z/L < 0), the following relations have been found:

φm = (1 − 16z/L)^−1/4 φ_h = (1 − 16z/L)^−1/2 and for the stable case (z/L > 0):

φm= φh= 1 + 5z L.

An analogous function for the water vapour profile is often set equal to the expression for φh.

(37)

Figure 4.5: A summary of dimensionless gradients of wind speed (a) and temperature (b) as they have been experimentally determined for the surface layer (From: Yaglom, blm 1977).

Upon integration of (4.1) and (4.2), we obtain

U (z) = u?

κ{ln(z/z₀) − Ψ_m(z/L)}

Θ(z) − Θ0 = θ?

κ{ln(z/z0) − Ψh(z/L)}

If z/L < 0, the expressions for Ψmand Ψhbecome

Ψm = 2 ln 1 + x 2

+ ln 1 + x² 2

− 2 tan⁻¹ x + π/2

Ψh = 2 ln 1 + y 2

, and otherwise (if z/L > 0)

Ψm= Ψh= −5z/L,

where x = (1 − 15z/L)^1/4 and y = (1 − 9z/L)^1/2(see figure 4.5).

If, for whatever practical reason, a power law for the wind speed profile is preferred, the exponent p is derived from the previous equations using (3.21):

p = φm(z12/L) ln z12/z0− Ψ(z12/L), where z12=√

z1z2.

Turbulence in the surface layer

(38)

The following (empirical) relations hold in a stable, flat surface layer:

σ_u/u_? = 2.39 ± 0.03 σv/u? = 1.92 ± 0.03 σ_w/u_? = 1.25 ± 0.03

In unstable conditions, this relation becomes

σw

u_? = 1.25

1 −3z

L

^1/3 .

In the unstable situation, the horizontal-velocity fluctuations are too strongly influenced by variations on a larger scale. The Monin-Obukhov theory starts to lose its validity. For the temperature variance we have

σθ/|θ?| = 0.95( z

−L)⁻¹³.

Expressions for the unstable situation, that also apply to the surface layer, will be given in section 4.4. (see also Garratt section 3.5)

Spectra

There exists a plethora of spectra and spectral data in the surface layer, for example in ?. For our purposes, it suffices to know that

f S_u(f )

u²_? = 102 n (1 + 33n)^5/3 f Sv(f )

u²_? = 17 n

(1 + 9.5n)^5/3 f Sw(f )

u²_? = 2.1 n (1 + 5.3n)^5/3,

(4.3) in which n is the dimensionless frequency (n = f z/U ).

4.3 The neutral boundary layer

The Ekman layer

In this section, we will consider a neutral, homogeneous and stationary boundary layer. Unlike previously, the Coriolis force will also be taken into account.

We will assume that the flow above the boundary layer is geostrophic, i.e. the

(39)

Figure 4.6: Force balance in the Ekman layer. P is the pressure gradient, Co the Coriolis force, and F r is friction.

pressure gradient and the Coriolis force balance each other. From the momentum equation, Eq. 3.3, neglecting the buoyancy term, it follows that

0 = f (V − Vg) −∂uw

∂z 0 = −f (U − Ug) −∂vw

∂z (4.4)

(see figure 4.6). Ug and Vg are the components of the geostrophic wind, by definition being equal to the wind speed while neglecting the turbulent stress in (4.4). These components are expressed as

U_g = − 1 ρ0f

∂Π

∂y

Vg = 1

ρ₀f

∂Π

∂x.

The Coriolis parameter, f , equals 2Ω sin φ, where φ represents the latitude, Ω the angular velocity (rad/s) of the Earth axis. In order to solve eqs. (4.4), K-theory is applied. Suppose that

uw = −Km

∂U

∂z

vw = −Km

∂V

∂z.

After substitution into (4.4), we obtain

0 = f (V − V_g) + K_m∂²U

∂z² 0 = −f (U − Ug) + Km

∂²V

∂z²,

BOUNDARY-LAYER METEOROLOGY