Klaus GierensInstitut für Physik der Atmosphäre

DLR Oberpfaffenhofen

Modelling of Cirrus Clouds

(MOD 10)

(MOD 11)

Overview

MOD11: Numerical modelling of important microphysical processes in cirrus clouds

MOD12: Stochastic cloud modelling

Klaus GierensInstitut für Physik der Atmosphäre

DLR Oberpfaffenhofen

Numerical modelling of important microphysical processes in cirrus clouds

MOD 11

Overview

• Problems special to cirrus modelling• Model types• Bulk microphysics models• Processes and their representation in my bulk model• Some modelling examples

Problems special to cirrus modelling - Radiation

1. Cirrus clouds may heat or cool the Earth-Atmosphere System depending onmicro-/macrophysical propertiestemperature (altitude)generation mechanism (incl. synoptic situation, geogr. location)

2. Complex ice crystal shapes (inter alia T- and Si- dependent) render calculation of radiative transfer a tough problem

Problems special to cirrus modelling – Ice formation

3. Various modes of ice crystal formationhomogeneous freezing of aqueous solution dropletsheterogeneous modes:

deposition freezing immersion freezing condensation freezingcontact nucleationand still other modes

from Vali, 2004

Problems special to cirrus modelling - Supersaturation

4. Cirrus clouds have only a loose relation to ice saturation, viz. they do not form at saturation once formed, they are not very strongly attracted by the

equilibrium state

Consequently:

there is plenty of ice supersaturated,

yet clear air in the UT

(sometimes marked by persistent contrails)

Cirrus clouds are embedded in

supersaturated air masses

RHi pdfs within cirrus have long tails into the

supersaturated regime Spichtinger et al., 2004

cloudy air

INCA data

Cirrus and Climate Change — an unsolved problem

Freezing/nucleation thresholds are high above saturation

extremal states in the RHi field

Extremal states react much more sensitive to changes of background conditions than do averages.

Hence it is difficult to estimate, how the probability will change in a changing climate that in the RHi field the nucleation thresholds will be surpassed.

Example:mean Si increases from 10 to 11%but probability to surpass 40%decreases by about 1/3wrt to the earlier pdf

Model types

• Models are a compromise between - numerical effort

o CPU costs, o computing and turnaround time, o memory and storage requirements

- and scientific ambition.

• Models with clouds usually combine - sophisticated dynamics with simple microphysics (NWP, GCM)

simple microphysics: bulk microphysics- simple dynamics with a elaborate microphysics; mostly box

models with size resolved microphysics: bin microphysics - trajectory calculations with single particle microphysics (recent

development)

Peculiar model types

• Models with both elaborated microphysics and detailed dynamics (e.g. Grabowski’s superparameterisation) are extremely expensive (in terms of computing power).

• Box models with bulk microphysics are almost never used. But they are very quick and one can learn a lot playing around with such a model (see Gierens, ACP, 2003).

Bulk microphysics models

• Bulk microphysics: balance equations only for few total concentrations that characterise a cloud. Typically- Mass concentration (1st moment of mass distribution)- Number density (0th moment)- Traditionally, many bulk models only transport the 1st moment

(e.g. the classical Kessler scheme)- Now, more dual-moment schemes (0th and 1st moments)

• Bulk schemes are usually used in NWP and GCM models and in many mesoscale models.

• Bulk schemes are much faster than bin microphysics schemes, at the expense of giving up information on size distribution (and probably also realism).

• The ECMWF model uses cloud coverage as a prognostic variable in addition to vapour and liquid/ice water concentrations.

Mathematical modelling of clouds, bulk version

• Needs assumption on probability density function type for the masses (or sizes) of the various hydrometeor and aerosol classes considered in the model.

• Marshall-Palmer (i.e. exponential) • gamma• log-normal • uni-modal, bi-modal, multi-modal

• Note: only the type of the pdf is chosen initially. The parameters generally change with time during evolution of the model cloud.

Common mass pdfs and their moments

avoid too many parameters!

• Number of parameters that fix the pdf should not exceed the number of prognostic variables by much.

• Parameters should be functions of the prognostic variables.

• The functional dependence should be understandable.

• It is difficult to determine a priori, how higher moments (skewness, curtosis, etc.) will evolve with the evolution of a cloud.

• Higher moments are difficult to determine from data- sensitive to outliers.

Processes to be included in a (pure) cirrus model

• Nucleation of the ice phase from- aerosol

o liquid (homogeneous)o solid (heterogeneous, various modes)

- water droplets• Crystal growth and evaporation• Crystal sedimentation• Crystal aggregation

• Aerosol dynamics and chemistry (parts of it implicit in nucleation)• Radiation (may feed back on growth/evap rates)

Processes marked in red are currently included in the Spichtinger/Gierens bulk cirrus physics of EuLag.

a typical set of equations

Equations used in the two-moment bulk cirrus scheme by Spichtinger and Gierens. Note also the two forms of ice!

Homogeneous nucleation of aqueous solution droplets

• Parameterisation after Koop et al.

• Critical supersaturation

• Nucleation rate J given as polynomial of awawi.

• In equilibrium the water activity equals the relative humidity wrt liquid water. Non-equilibrium occurs in strong updraughts.

• Integration over droplet size distribution:

- Actual droplet volume derived by inversion of Köhler equation.

- For a log-normally distributed dry aerosol mass Gauss-Hermite integration works fine (Gierens and Ström, JAS, 1998).

Köhler equation

For a given ambient relative humidity the equilibrium size of a solution droplet is given by the Köhler equation. Simplest form:S = A/r B/r3

A/r is the Kelvin term,B/r3 is the Raoult term

Heterogeneous nucleation

• simplest assumption possible:- a certain number of solid aerosol particles (typically 1 to 50 per

cubic centimetre)- freeze to ice at a certain supersaturation (typically 130%).

• On evaporation of het. ice, these aerosols are set free and can form new ice afterwards.

Deposition growth and evaporation

• Parameterisation after Koenig (JAS, 1971):- dm/dt = a mb with temperature, pressure, and

supersaturation dependent coefficients a,b. • Corrections for kinetic growth regime (small ice crystals)

Integration over mass distribution

• Integration over mass distribution:- Diffusion regime:

- dIWC/dt = a µb

- Kinetic regime:

- dIWC/dt = a µb+ / m0

• a > 0 implies growth of the ice mass concentration: dIWC/dt > 0- the ice number density is the constant.

• a<0 implies crystal evaporation (dIWC/dt < 0).- the ice number concentration decreases then, but with a higher

relative rate than the ice mass:

- (Nt-1 Nt)/ Nt-1 = [(IWCt-1 IWCt)/ IWCt-1] with =1.1

- (Harrington et al., 1985)

Two-moment sedimentation scheme

• Flux densities for ice mass and number concentrations

• Empirical relation between crystal mass and terminal velocity

Two-moment sedimentation scheme, cont’d

• Allows to express mass and number related terminal velocities as:

• Since large crystals fall faster than small ones, one must have

- vt,m > vt,n

- in other words: µ+1µ0 > µ µ1.

• This inequality is always fulfilled (Gierens and Spichtinger, SPL, subm.)

Simulation of different sedimentation

• Shape of ice crystals: columns

• Initialising of a thin cirrus cloud at t=0s ( IWC = 10 mg / m3 , N = 100 / dm3 ) in the altitude range 8.5 -9.5 km

• Simulation time: t = 3600 s

viwc vnc viwc = vnc

The two-moment sedimentation scheme nicely obyes the principle that large crystals fall faster than smaller ones. This is not so in the one-moment scheme.

Two vs. One-moment sedimentation, other effects

Effect on vertical distribution ofice water mass and numberconcentration.

Effect on SW and LWextinction per modellayer. Larger optical thickness in the 2-moment scheme.

Sensitivity studies: homogeneous vs. heterogeneous nucleation

Results from the DFG project „Dünner Zirrus“ (thin cirrus).

Setup for idealised 2D simulations

Model domain:

horizontal resolution dx = 100 m, horizontal extension: 6.3 km

vertical resolution dz = 50 m, vertical extension: 6 km, i.e. 6-12 km

time step dt = 1s, simulation time 6 h = 21600 s

constant vertical motion for whole model domain (i.e. adiabatic cooling)

w = 3 / 4.5 / 6 cm/s

• Set of number densities of ice nuclei: Ni = 1 / 3 / 5 / 7 / 10 / 30 / 50 L-1

• Set of thresholds for heterogeneous nucleation: RHihet = 110 / 130 / 140 %

• Additional temperature fluctuations: T = 0.1 / 0.05 / 0.01 / 0.005 / 0.001 K

Start profiles

Varying ice nuclei number density Ni

In the following mean values over all 64 columns are shown:

• x-axis: time in minutes

• z-axis: altitude in metres

• colour bar: relative humidity with respect to ice

Isolines of equal ice crystal number densities

purple: ice crystals formed by homogeneous nucleation

black: ice crystals formed by heterogeneous nucleation

Ni = 1L-1, w = 4.5 cm/s, RHihet = 130 %

Time (min)

Alti

tude

(m

)

Ni = 3L-1, w = 4.5 cm/s, RHihet = 130 %

Alti

tude

(m

)

Time (min)

Ni = 5L-1, w = 4.5 cm/s, RHihet = 130 %

Alti

tude

(m

)

Time (min)

Ni = 7L-1, w = 4.5 cm/s, RHihet = 130 %

Alti

tude

(m

)

Time (min)

Ni = 10L-1, w = 4.5 cm/s, RHihet = 130 %

Alti

tude

(m

)

Time (min)

Ni = 30L-1, w = 4.5 cm/s, RHihet = 130 %

Alti

tude

(m

)

Time (min)

Ni = 50L-1, w = 4.5 cm/s, RHihet = 130 %

Alti

tude

(m

)

Time (min)

Results 1

• If one of these competing nucleation mechanisms (heterogeneous/homogeneous) can produce many ice crystals, relative humidity can be reduced effectively.

• Two different regimes:- few heterogeneous ice nuclei: homogeneous nucleation is effective- many heterogeneous ice nuclei: heterogeneous nucleation is

effective• between these two regimes the cloud is very sensitive to the number of

ice nuclei; often there is persistent ice supersaturation within the simulated clouds, reaching rather high values.

• transition between the two regimes depends on the relation between three time scales:

growth - sedimentation - cooling

varying threshold humidity for heterogeneous nucleation

• In the regimes where one formation mechanism is dominant only marginal changes are due to different thresholds

• In the range where no process is dominant a change in the threshold affects the properties of the clouds quite seriously

• For low thresholds a “secondary cloud formation” is observed:- Ice crystals sediment and evaporate in the sub saturated layers

below the cloud

Moistening of the sub saturated layer

Collection of aerosols in this layer

Due to cooling cloud formation by heterogeneous nucleation

Ni = 5L-1, w = 4.5 cm/s, RHihet = 110 %

Alti

tude

(m

)

Time (min)

Secondary cloud formation

Ni = 5L-1, w = 4.5 cm/s, RHihet = 130 %

Alti

tude

(m

)

Time (min)

Ni = 5L-1, w = 4.5 cm/s, RHihet = 140 %

Alti

tude

(m

)

Time (min)

Varying ice nuclei number density Ni with additional temperature fluctuations

In the following mean values over all 64 columns are shown:

• x-axis: time in minutes

• z-axis: altitude in metres

• colour bar: relative humidity with respect to ice

Isolines of equal ice crystal number densities

purple: ice crystals formed by homogeneous nucleation

black: ice crystals formed by heterogeneous nucleation

Temperature fluctuations: Gaussian, T = 0.1 K

Ni = 1L-1, w = 4.5 cm/s, RHihet = 130 %, T=0K

Ni = 1L-1, w = 4.5 cm/s, RHihet = 130 %, T=0.1K

Ni = 7L-1, w = 4.5 cm/s, RHihet = 130 %, T=0K

Ni = 7L-1, w = 4.5 cm/s, RHihet = 130 %, T=0.1K

Ni = 50L-1, w = 4.5 cm/s, RHihet = 130 %, T=0K

Ni = 50L-1, w = 4.5 cm/s, RHihet = 130 %, T=0.1K

Varying IN number density Ni with temperature fluctuations

• In the regimes where one formation mechanism is dominant only marginal changes are due to temperature fluctuations

• In the range where no process is dominant temperature fluctuations affect the properties of the clouds quite seriously

• The effect is in two directions: Temperature fluctuations can

o enforce the reduction of relative humidityo slow down the reduction of relative humidity

Klaus GierensInstitut für Physik der Atmosphäre

DLR Oberpfaffenhofen

Stochastic cloud modelling

MOD 12

Stochastic cloud modelling (statistical schemes)

Main problem here: Parameterisation of cloud fraction (i.e. fractional cloud cover).

Problem for large scale models, not for cloud resolving models.

In a CRM a grid box is either cloudy or cloud free (binary or 0-1 scheme). Some old GCMs also use this binary assumption of total or zero cloud cover.

The 0-1 schemes neglect sub-grid variability. This leads to errors in all computations, where quantities depend nonlinearly on liquid or ice water path or concentration.

Statistical cloud schemes would allow to consistently treat sub-grid variability in cloud microphysical processes and in radiation.

A cloud resolving model runs reasonably well with a 0-1 scheme

In a large scale model the results of a 0-1 scheme are unsatisfying

A cloud fraction looks somewhat better, although problems of cloud overlap assumptions arise, in particular for radiative transfer. Most (all?) models do not assume a variable cloud fraction in the vertical within one grid layer.

RH-controlled parameterisation of cloud cover

In some GCM schemes cloud cover is parameterised as a function of relative humidity, e.g. the so-called Sundqvist scheme of ECHAM.

Simple statistical scheme:

Clouds already form at Uc>100%, i.e. at sub-saturated conditions.

Interpretation:

Fluctuations of RH in the grid box Supersaturation somewhere clouds form in a fraction of the box.

Some schemes use also vertical wind speed to parameterise C.

21

1

),max(11

c

c

U

UUC

Statistical schemes

Working principle

Consider a phase space (T,RH). In a certain part of the phase space clouds can form, in the remaining part not.

Examples:

Water clouds: RH>100%

Ice clouds: RHi > RHicrit (T)

The model predicts at every time step and for each grid box a mean state <(T,RH)>.

Statistical schemes, cont’d

If we know the probability density function of fluctuations of the phase point around the grid-box mean value, we can compute, how probable it is that a fluctuation reaches into the supercritical regime.

I call this probability the Overlap Integral .

can be interpreted as the actual cloud coverage C.

For numerical reasons it might be better to compute d / dt and from that dC/dt.

Examples for homogeneous nucleation and contrail formation

Red line: critical supersaturation for homogeneous nucleation(Koop theory).Green dots: fluctuations of temperature and relative humidity around the grid mean state(-50°C, 140%).

= (number of dots above the red line) / (total number of dots)

Phase diagram for formation of persistent contrails(for two pressure levels).

Probability density functions for fluctuations

Problem:• how are fluctuations of the phase state variables distributed.

• The distribution may depend on - location, - time (e.g. season), - in particular on the spatial scale

o spatial resolution of the modelo spatial resolution of data (correlation lengths)

Probability density functions for fluctuations, cont’d

• Generally, pdfs are chosen in an ad hoc way, - data on fluctuations almost non-existent. - pdf selection according to criteria outside of physics, - more inside of mathematics and numerics.- symmetrical pdfs often used, BUT- symmetric pdfs cannot be the true nature of the fluctuations

since temperature and relative humidity (or other humidity variables) cannot be negative.

• Apart from measurements, distributions of fluctuations are also sometimes obtained from cloud resolving model runs e.g. Adrian Tompkins). It is clear that many runs are needed to get a good statistical ensemble.

measured statistics of instantaneous fluctuations

MOZAIC data (one year)Gierens et al., Ann. Geophys., 1997

Analytical formulation of the fluctuations

The measured fluctuations (on a T42 grid scale, i.e. 250×250 km2) follow closely a Cauchy distribution (Lorentz line shape):

() = ( / ) / (2 + 2)

Cauchy distribution: no moments! (not even a mean value). Widely extended tails.

BUT:convolution of two Cauchy distributions yield another Cauchy distribution

1() 2() = 1+2()

joint probability density of (T,RH)-fluctuations

linearise RH(T)

sum of two random variables convolution

evaluate convolution integral

joint probability density of (T,RH)-fluctuations, cont’d

Insert the two Cauchy distributions ().

Result is a product of the original Cauchy distribution for Twith a “rotated” Cauchy distribution for RH + AT

Theoretical joint pdf of (T,RH) fluctuations

joint pdf of (T,RH) fluctuations constructed from MOZAIC data

Overlap integral for contrail formation, analytically

The calculation of theoverlap integral effectivelysmears out the boundariesin the phase space.

overlap integral for Koop parameterisation, numerically

practical considerations – random numbers with certain pdf

• Random number generators usually produce uniform distribution of variable R on [0,1).

• Cauchy distribution: distribution of tan , with uniformly distributed within [/2, /2].

• hence set random = tan (R - /2)

• For other distributions: inversion of cumulative distribution function F (integral of the pdf)

R

F

0

1

x

random number x is:

x = F-1 (R)

where R is a randomnumber in [0,1) producedby a generator.

practical considerations - d/dt

• total derivative- d/dt = (/T) (dT/dt) + (/RH) (dRH/dt)

• How to compute the partial derivatives of wrt the phase variables?

• Analytical expression (at least with 2D-Cauchy distribution) are VERY complex, unfeasible…

• Numerical approximation.• (/T) [(T) (T+dT)] / dT yields noisy results and needs a

lot of computing time for computing the random numbers.

• Better idea: see next slide!

practical considerations - d/dt, cont’d

Temperature derivative:shift the red line in ±T-direction by ±dT, count the number of points between the black lines, divide by total number of points and divide by 2 dT.

RH derivative:shift the red line in ±RH-direction by ±dRH, count the number of points between the black lines, divide by total number of points and divide by 2 dRH.