Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Investigation into nonlocal turbulence-closure at higher statistical order Modzelewski, Henryk 2004

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_2004-931579.pdf [ 13.5MB ]
Metadata
JSON: 831-1.0052396.json
JSON-LD: 831-1.0052396-ld.json
RDF/XML (Pretty): 831-1.0052396-rdf.xml
RDF/JSON: 831-1.0052396-rdf.json
Turtle: 831-1.0052396-turtle.txt
N-Triples: 831-1.0052396-rdf-ntriples.txt
Original Record: 831-1.0052396-source.json
Full Text
831-1.0052396-fulltext.txt
Citation
831-1.0052396.ris

Full Text

I N V E S T I G A T I O N INTO N O N L O C A L T U R B U L E N C E - C L O S U R E A T H I G H E R S T A T I S T I C A L O R D E R by Henryk Modzelewski M . S c , The Warsaw University, Poland, 1989 A THESIS S U B M I T T E D IN P A R T I A L F U L L F I L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F D O C T O R O F P H I L O S O P H Y in T H E F A C U L T Y O F G R A D U A T E STUDIES Department of Earth and Ocean Sciences Faculty of Science Atmospheric Science Programme We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F B R I T I S H C O L U M B I A July 2004 ©Henryk Modzelewski, 2004 Abstract This dissertation investigates whether higher-statistical-order nonlocal turbulence-closure can be utilized to describe some of the complex features of atmospheric turbulent flows. Transilient turbulence theory, a nonlocal forcing and mixing model, was chosen as a framework. This theory describes turbulent transport as a two-dimensional 'transilient' matrix of exchange coefficients between different layers in the fluid. A new higher-statistical-order nonlocal parameterization of the transilient matrix is proposed. The studies of former parameterizations of transilient matrices have indicated the following main deficiencies: lack of asymmetry of upward and downward motions in convective boundary layers, improper vertical distribution of turbulent mixing, ad hoc definition of diagonal elements of transilient matrices, and unrealistically large ranges of turbulent transport. This dissertation proposes a parameterization that addresses these deficiencies. The following concepts are utilized in the new parameterization: a virtual turbulent-transport eddy, convective potential-energy and shear potential-energy of mean-flow instabilities, decom-position of forcing into symmetric and asymmetric components, prognostic equations for the nonlocal turbulence kinetic energy budget and for the 3rd statistical moment of the vertical ve-locity, a limited turbulent-transport range, and a transilient tendency-matrix as a precursor of the transilient matrix. The parameterization utilizes eight parameters. A testbed atmospheric boundary-layer model is created to investigate, evaluate, and compare the new closure with three other turbulence-closure models and with two large-eddy-simulation data sets. The evaluation of testbed model led to the following conclusions. The parameterization of asymmetry qualitatively predicts the overall asymmetry of turbulent transport. A n improper ver-tical distribution of mixing in convective boundary layers is evident, caused by improper height and scale-wise structure of buoyancy forcing. The utilization of the new transilient tendency-matrix successfully removed the necessity to parameterize the diagonal elements of transilient matrices. The limited transport-range approximation causes sensitivity discretization. Finally, the parameterization is computationally expensive, mainly due to the cost of solving the prognos-tic equations. In summary, the new parameterization has demonstrated the potential to improve turbulence parameterization, but more work needs to be done. ii Table of Contents Abstract ii Table of Contents iii List of Tables vi List of Figures vii Notation, symbols, and abbreviations xi Preface xviii Acknowledgments xix 1 Review of turbulence properties 1 1.1 Properties of turbulent flows 1 1.1.1 Prediction of turbulence 2 1.1.2 Spectral properties of turbulence 6 1.1.3 Properties of energy-containing eddies 7 1.2 Properties of turbulent transport 8 1.2.1 General properties 8 1.2.2 Specific features of convective transport 9 2 Review of turbulence modeling 11 2.1 The turbulence closure problem 11 2.2 Review of turbulence closures 11 2.2.1 Overview 11 2.2.2 Local closures 13 2.2.3 Nonlocal closures 14 3 Review of Transilient Turbulence 17 3.1 Fundamentals of Transilient Turbulence Theory 18 3.1.1 Transilient matrix concept . 19 3.1.2 Turbulence statistics based on the transilient matrix 24 3.2 Parameterizations of the Transilient Matrix 26 3.2.1 Transilient matrix based on Richardson number 27 3.2.2 Transilient matrix based on diagnostic T K E 28 3.3 Summary 32 3.3.1 The motivation for this research 33 ii i 4 A prognostic model for the transilient matrix 35 4.1 Prognostic transilient model (PTM) concepts 36 4.1.1 Turbulent-transport mean-flow state operators 36 4.1.2 Virtual turbulent-transport eddy concept 37 4.1.3 Extension of nonlocal static stability criteria 39 4.1.4 Conversion between potential and kinetic energy 42 4.2 Prognostic equations for W2 and W3 43 4.2.1 Turbulence generation by mean-flow instabilities 45 4.2.2 Turbulence energy-cascade 50 4.2.3 Asymmetry decay 52 4.2.4 Summary 52 4.3 Mixing potential 53 4.3.1 A n M P parameterization 53 4.3.2 Mass-flux conservation algorithm 54 4.4 Transilient matrix 55 4.4.1 Transilient tendency-matrix 55 4.4.2 Transforming [C](t) into [C](t,At) 56 4.5 Notes on T3 statistics 56 4.6 Summary 56 5 Estimation of P T M parameters 59 5.1 Available L E S cases 60 5.2 P T M experiment design 61 5.3 Parameter-space exploration procedure 62 5.4 Sensitivity case-studies 63 5.4.1 E K M case - idealized Ekman layer 63 5.4.2 00WC case - neutral boundary layer 66 5.4.3 24F case - convective boundary layer 68 6 P T M evaluation and comparison 72 6.1 P T M comparison with other models for the LES-database cases 72 6.1.1 Neutral layer with strong capping inversion - 00SC case 72 6.1.2 Weak capping inversion - 05WC case 75 6.1.3 Strong capping inversion - 03SC, 05SC, and 24SC cases 77 6.1.4 Baroclinic atmosphere - 15B and 24B cases 82 6.1.5 Summary 85 6.2 Ebert case 86 6.3 Sensitivity to discretization 92 6.3.1 Sensitivity to integration time-step 93 6.3.2 Sensitivity to vertical resolution 96 6.3.3 Summary 100 6.4 Conclusions 100 7 Conclusions and recommendations for future work 102 References 105 A Nonlocal closures 111 iv B L E S simulations 114 B . l L E S database for P B L model evaluation 114 B . l . l L E S model 114 B.1.2 L E S database at N C A R 115 B . l . 3 Creation of the L E S database 116 B.2 Ebert case 117 C P B L models used in comparison study 119 D Numerical procedures and approximations 121 D . l P T M equations step-by-step 121 D.2 Numerical approximations 122 v List of Tables 2.1 Summary of turbulence closures (after Stull (1994b)). 0 denotes mean potential temperature, 9 denotes perturbation of mean potential temperature, w denotes perturbation of mean vertical velocity, dt denotes time derivative, and dz denotes height derivative 12 5.1 Archetypical P B L cases of L E S runs, the output from which are used to compare with P T M and with other P B L models. (The barotropic cases have imposed constant geostrophic wind with height, while the baroclinic cases impose shear in the geostrophic wind. A l l cases are barotropic, unless noted otherwise.) 60 6.1 Combined statistics for all compared models. Smaller values of C R M S and A C R M S , mean less difference from the LES simulations. Higher values of C S R C and A C S R C indicate better correlation with L E S . Values of C M E and A C M E are meaningful only when C R M S and A C R M S are comparable, and then smaller absolute values indicate less difference from the L E S results 73 6.2 Turbulence invariant, Im, and asymmetry indices, A\ through A4, for transilient matrices predicted by P T M and measured by Ebert et al. (1989). The indices with values smaller than 1 0 - 1 0 were replaced with zeros 91 B . l External forcing parameters for L E S simulations and P B L runs: Ug and Vg are the components of geostrophic wind, ZQ is surface roughness-length, and wd is the surface sensible heat-flux. dzQ is the initial gradient of potential temperature in the free atmosphere and CI indicates if strong capping inversion was imposed in the initial conditions of LES simulation. (Ayotte et al., 1996) 115 B.2 Internal scales of LES simulations (large-eddy overturning time, r , and boundary-layer depth, Zi) used to determine vertical resolution and integration time-step of P B L runs. Also given is the length, tSim, of P B L runs. (Ayotte et al., 1996) 116 vi List of Figures 1.1 A n idealized profile of mean virtual potential temperature, 0^, and vectors showing magnitude and direction of turbulent heat-flux, w9v, in the convective boundary layer (CBL) . z is height above ground. The turbulent C B L consists of three distinct sub-layers: a super-adiabatic layer in the lower part, an adiabatic layer in the middle, and a sub-adiabatic layer in the upper part. Positive turbulent heat-flux typically extends into the upper part of the C B L 5 3.1 Vertical cross-section of a single column of a typical medium-resolution mesoscale-forecast-model grid. The thick-line boxes represent scopes of circulation of typical convective eddies (with horizontal to vertical aspect ratio of roughly 2) with various velocity scales as indicated by w, preserving approximately the aspect ratio of grid sizes 19 3.2 Graphical representation of the transilient coefficients for a 5-layer model. The coefficients define the relative contribution of eddies of different sizes to the total turbulent transport between different layers of the model 20 3.3 Example of a transilient matrix for a 5-layer model with two turbulent sub-domains separated by a non-turbulent layer. ' D ' indicates destination index and 'S' indicates source index. A n arrow in the lower-right corner indicates the main diagonal of the matrix. The turbulent subdomains are indicated by thick-line boxes. A non-turbulent domain is represented by an oval-finished cross. The matrix is flipped, as compared to the typical mathematical notation, so that the main diagonal has the same direction, but the lowest indices are in lower-right corner. The notation was introduced by Stull (2000) for easier interpretation of a Transilient Matrix; namely, higher in the matrix corresponds to higher in the atmosphere 21 3.4 Limits of the transilient matrix, [C](At) from Figure 3.3, as At approaches 0 or oo. 23 4.1 A n illustration of nonlocal stability detection using Stuff's vs. the generalized-C A P E algorithm. A warmer layer is introduced into sub-adiabatic virtual potential temperature profile. Vertical lines represent the depth of the convection-driven turbulent layer as detected using: Stuff's nonlocal static stability method (limited by stars); and the generalized-CAPE method (limited by boxes). These are the statically unstable regions 41 vn 5.1 P T M m a t r i c e s for the E K M case: [TVS], [Y], a n d [C] , w h e r e s is t h e h e i g h t o f a sou rce layer , a n d d is t h e he igh t o f a d e s t i n a t i o n l ayer . T h e f o l l o w i n g s h a d -i n g / c o n t o u r i n g schemes are u s e d t o d i s p l a y t h e m a t r i c e s for a l l L E S - d a t a b a s e cases. T e n e q u a l l y s p a c e d s h a d i n g levels are u s e d for [7\S] a n d [AS] w i t h a d e l i m i t e r ev-e r y 0.1 (i .e. 0 . 1 , . . . ,0.9) o f difference b e t w e e n t h e m a x i m u m a n d m i n i m u m va lues . E l e v e n s h a d i n g levels are used for [Y] a n d [C], a n d d e l i m i t e r s are set t o l / 2 n , w i t h n = 1 , . . . , 10, o f the m a x i m u m v a l u e (i .e. , 0 .5, 0 .25, 0 .125, 0 .0625, 0 .03125, 0 .015625, 0 .0078125, 0 .00390625, 0 .00195313, a n d 0 .000976563 o f t h e m a x i m u m va lue ) i n o r d e r t o s h o w be t t e r t h e i r p r o p e r t i e s for s h o r t i n t e g r a t i o n t i m e - s t e p s . T h e s h a d i n g d a r k n e s s increases w i t h t h e va lue , f r o m n o s h a d i n g for s m a l l va lues t o b l a c k for l a rge ones 64 5.2 V e r t i c a l p rof i les for t h e E K M case: P T M ( s o l i d l i ne ) v s . L E S ( d a s h e d l ine ) 65 5.3 P T M m a t r i c e s for t h e 0 0 W C case: [DS], [Y], a n d [C]. See c a p t i o n i n F i g u r e 5.1 for p l o t t i n g de t a i l s 66 5.4 V e r t i c a l prof i les for t h e 0 0 W C case: P T M (so l id l i ne ) v s . L E S ( d a s h e d l i n e ) , [au] s t a n d s for a r b i t r a r y u n i t s 67 5.5 P T M m a t r i c e s for t h e 2 4 F case: [7\S], [AS], [Y], a n d [C] . See c a p t i o n i n F i g u r e 5.1 for p l o t t i n g de t a i l s 68 5.6 V e r t i c a l prof i les for t h e 2 4 F case: P T M ( s o l i d l i ne ) v s . L E S ( d a s h e d l i n e ) , [au] s t a n d s for a r b i t r a r y u n i t s 69 5.7 2 4 F case: P T M prof i les for t h e s y m m e t r i c s i m u l a t i o n ( l o n g - d a s h l i n e ) , a n d a s y m -m e t r i c s i m u l a t i o n ( s h o r t - d a s h l ine) v s . L E S ( s o l i d l i n e ) , [au] s t a n d s for a r b i t r a r y u n i t s . A l s o s h o w n are [715]* w i t h s y m m e t r i c f o r c i n g on ly , s i de -by - s ide w i t h [T\S] i n c l u d i n g a s y m m e t r i c f o r c i n g . See c a p t i o n i n F i g u r e 5.1 for p l o t t i n g d e t a i l s 70 6.1 V e r t i c a l p rof i les for t h e 0 0 S C case: P T M ( s o l i d l ine) v s . L E S ( d o t t e d l i n e ) , T R A N ( s i n g l e - d o t t e d d a s h e d l i n e ) , K P R O ( d o u b l e - d o t t e d d a s h e d l i n e ) , a n d M Y 2 5 ( t r i p l e -d o t t e d d a s h e d l i n e ) , [au] s t a n d s for a r b i t r a r y u n i t s 74 6.2 P T M m a t r i c e s for t h e 0 0 S C case: J\S, Y , a n d C, w h e r e s is t h e h e i g h t o f a sou rce l ayer , a n d d is t h e he igh t o f a d e s t i n a t i o n l ayer . T h e f o l l o w i n g s h a d i n g / c o n t o u r i n g schemes are u s e d t o d i s p l a y t h e m a t r i c e s for a l l L E S - d a t a b a s e cases. T e n e q u a l l y s p a c e d s h a d i n g levels are u s e d for [7\S] a n d [AS] w i t h a d e l i m i t e r e v e r y 0.1 (i .e. 0 . 1 , . . . ,0.9) o f dif ference b e t w e e n t h e m a x i m u m a n d m i n i m u m va lues . E l e v e n s h a d -i n g leve ls are u s e d for [Y] a n d [C], a n d d e l i m i t e r s are set t o l / 2 n , w i t h n = 1 , . . . , 10, o f t h e m a x i m u m v a l u e (i .e. , 0.5, 0 .25, 0.125, 0 .0625, 0 .03125, 0 .015625, 0 .0078125, 0 .00390625, 0 .00195313, a n d 0 .000976563 o f t h e m a x i m u m va lue ) i n o r d e r t o s h o w b e t t e r t h e i r p r o p e r t i e s for s h o r t i n t e g r a t i o n t ime- s t eps . T h e s h a d i n g d a r k n e s s i n -creases w i t h t h e va lue , f r o m no s h a d i n g for s m a l l va lues t o b l a c k for l a rge ones . . . 75 6.3 V e r t i c a l prof i les for t h e 0 5 W C case: P T M ( s o l i d l ine) vs . L E S ( d o t t e d l i n e ) , T R A N ( s i n g l e - d o t t e d d a s h e d l i n e ) , K P R O ( d o u b l e - d o t t e d d a s h e d l i n e ) , a n d M Y 2 5 ( t r i p l e -d o t t e d d a s h e d l i n e ) , [au] s t a n d s for a r b i t r a r y u n i t s 76 6.4 P T M m a t r i c e s for the 0 5 W C case: 1\S, AS, Y , a n d C . See c a p t i o n i n F i g u r e 6.2 for p l o t t i n g de t a i l s 77 6.5 V e r t i c a l p rof i les for t h e 0 3 S C case: P T M ( so l id l ine) v s . L E S ( d o t t e d l i n e ) , T R A N ( s i n g l e - d o t t e d d a s h e d l i n e ) , K P R O ( d o u b l e - d o t t e d d a s h e d l i n e ) , a n d M Y 2 5 ( t r i p l e -d o t t e d d a s h e d l i n e ) , [au] s t a n d s for a r b i t r a r y u n i t s 78 6.6 V e r t i c a l p ro f i l e s for t h e 0 5 S C case: P T M ( s o l i d l ine) vs . L E S ( d o t t e d l i n e ) , T R A N ( s i n g l e - d o t t e d d a s h e d l i n e ) , K P R O ( d o u b l e - d o t t e d d a s h e d l i n e ) , a n d M Y 2 5 ( t r i p l e -d o t t e d d a s h e d l i n e ) , [au] s t a n d s for a r b i t r a r y u n i t s 79 v i i i 6.7 Vertical profiles for the 24SC case: P T M (solid line) vs. L E S (dotted line), T R A N (single-dotted dashed line), K P R O (double-dotted dashed line), and MY25 (triple-dotted dashed line), [au] stands for arbitrary units 80 6.8 P T M matrices for the 03SC case: TVS, AS, Y , and C. See caption in Figure 6.2 for plotting details 81 6.9 P T M matrices for the 05SC case: TiS, AS, Y , and C. The grey background in [AS] corresponds to value of zero. See caption in Figure 6.2 for plotting details. . 81 6.10 P T M matrices for the 24SC case: TiS, AS, Y , and C. The grey background in [AS] corresponds to value of zero. See caption in Figure 6.2 for plotting details. . 81 6.11 Vertical profiles for the 15B case: P T M (solid line) vs. L E S (dotted line), T R A N (single-dotted dashed line), K P R O (double-dotted dashed line), and MY25 (triple-dotted dashed line), [au] stands for arbitrary units 83 6.12 Vertical profiles for the 24B case: P T M (solid line) vs. L E S (dotted line), T R A N (single-dotted dashed line), K P R O (double-dotted dashed line), and MY25 (triple-dotted dashed line), [au] stands for arbitrary units 84 6.13 P T M matrices for the 15B case: TVS, AS, Y , and C . The grey background in [AS] corresponds to value of zero. See caption in Figure 6.2 for plotting details. . 85 6.14 P T M matrices for the 24B case: TVS, AS, Y , and C. The grey background in [AS] corresponds to value of zero. See caption in Figure 6.2 for plotting details. . 85 6.15 P T M matrices for the Ebert case: [TVS] and [AS], where s is the height of a source layer, and d is the height of a destination layer. Twenty equally spaced shading levels are used with a delimiter every 0.05 (i.e., 0.05,.. .,0.95) of difference between the maximum and minimum values. The shading darkness increases with the value, from no shading for small values to black for large ones. The grey background in [AS] corresponds to value of zero 86 6.16 P T M vs. Ebert [C] m for te = 0.02r, where s is the height of a source layer, and d is the height of a destination layer. Twenty equally spaced shading levels are used with a delimiter every 0.05 (i.e., 0.05,.. .,0.95) of difference between the maximum and minimum values. The shading darkness increases with the value, from no shading for small values to black for large ones. The grey background in Ebert [C] m corresponds to value of zero, and the small negative values of in Ebert's matrix (plotted as brighter than background grey) are artifacts of the numerical accuracy of the L E S model 87 6.17 P T M vs. Ebert [C] m for te = 0.2r. The grey background in Ebert [ C ] m corresponds to value of zero, and the small negative values of in Ebert's matrix (plotted as brighter than background grey) are artifacts of the numerical accuracy of the L E S model. See caption in Figure 6.16 for plotting details 88 6.18 P T M vs. Ebert [C] m for te = 0.5r. See caption in Figure 6.16 for plotting details. . 88 6.19 P T M vs. Ebert [C]m for te = 0.6r. See caption in Figure 6.16 for plotting details. . 89 6.20 P T M vs. Ebert [C] m for te = l . r . See caption in Figure 6.16 for plotting details. . 89 6.21 P T M vs. Ebert [C] m for te = 2.r. See caption in Figure 6.16 for plotting details. . 90 6.22 P T M vs. Ebert [C]m for te = 4.r. See caption in Figure 6.16 for plotting details. . 90 6.23 E K M case: L E S profiles (solid line), and P T M profiles for At equal T / 2 0 (short-dash line), T / 1 0 (medium-dash line), and r / 5 (long-dash line) 93 6.24 00WC case: L E S profiles (solid line), and P T M profiles for At equal r /20 (short-dash line), r /10 (medium-dash line), and r / 5 (long-dash line), [au] stands for arbitrary units 94 ix 6.25 24F case: L E S profiles (solid line), and P T M profiles for A t equal r /20 (short-dash line), r /10 (medium-dash line), and r /5 (long-dash line), [au] stands for arbitrary units 95 6.26 E K M case: L E S profiles (solid line), and P T M profiles for Az equal Zi/20 (short-dash line), Zi/10 (medium-dash line), and Zi/5 (long-dash line) 97 6.27 OOWC case: LES profiles (solid line), and P T M profiles for Az equal Zi/20 (short-dash line), Zi/10 (medium-dash line), and Zj/5 (long-dash line), [au] stands for arbitrary units 98 6.28 24F case: L E S profiles (solid line), and P T M profiles for Az equal zt/20 (short-dash line), Zi/10 (medium-dash line), and Zi/5 (long-dash line), [au] stands for arbitrary units 99 x Notation, symbols, and abbreviations Tensor notation The following common tensor notation is used throughout this work. Discrete functions and variables, i.e. matrices and vectors, are indicated with square brackets to distinguish them from continuous variables. To distinguish tensor indices from other subscripts, the i-th element of vector [V] is defined as V^y such that the i-th row and j-th column element of matrix [M] is M[j j]. For example [V] represents L V[N] and [M] represents [M] = Mi [Li] M< l,N] Mr, L[N,l] • • • M{N,N] J The arguments of the function are listed in parenthesis, e.g. [V](i) represents a vector dependent on time, while Vyj (t) represents j-th element of a vector dependent on time. Derivatives A shortened notation is used for partial derivatives; e.g., dz = Symbols and abbreviations Symbols and abbreviations are given in the following order: 1. lowercase Greek letters 2. uppercase Greek letters 3. Roman letters in alphabetical order (a) lowercase letters (b) uppercase letters (c) calligraphic letters (d) double-struck letters List of symbols and abbreviations: xi Lowercase Greek B parameter in off-diagonal elements of Y 7c DeardorfT's correction coefficient 7s convective counter-gradient correction 6ij Kronecker's 5 symbol S(z, Z) Dirac's 8 distribution e dissipation rate 6 fluctuation of mean potential temperature 0 0V fluctuation of mean virtual potential-temperature ©„ K wavenumber associated with scale I KM wavenumber of the largest turbulent eddies v kinematic viscosity of the fluid £ fluctuation of generic mean conserved mean-quantity S Q mass density T large-eddy overturning time Upper-case Greek A finite-difference operator Aa; horizontal size of a grid cell A i integration time-step of discrete model Az depth of a layer in finite-difference model A'° nonlocal difference operator 0 mean potential temperature 0 „ mean virtual potential-temperature Qy virtual potential temperature in the cloud 0 ^ virtual potential temperature in the surrounding environment @ v vertical mean of virtual potential-temperature A length-scale parameter A l J nonlocal average operator H generic conserved mean-quantity Es mean-quantity scaled into interval (0,1) \P passive tracer A au arbitrary units a" coefficient in linear mapping into (0,1) A horizontal area of V T T E Ai 1-st index of absolute spectral asymmetry A2 2-nd index of absolute spectral asymmetry A3 3-rd index of spectral asymmetry A4 4-th index of spectral asymmetry AS S velocity scale B b fluctuation of passive tracer B 6" coefficient in linear mapping into (0,1) B mean passive tracer B volume of V T T E Bdown downdraft sub-volume of V T T E xii Bup updraft sub-volume of V T T E Bo dimensionless constant B R N bulk Richardson number for convective-systems BRN bulk Richardson number for convective-systems C c fluctuation of passive tracer C C mean passive tracer Cp specific heat of air C A P E convective available potential-energy CAPE convective available potential-energy CAPE* generalized convective available potential-energy C B L convective boundary layer C F L Courant-Friedrichs-Levy condition for computational stability C P E convective potential-energy CPE convective potential-energy C transilient matrix <Cm mass-normalized transilient matrix subdomain transilient-matrix C transilient tendency-matrix Cd intermediate matrix in derivation of C Ck intermediate matrix in derivation of C D d destination-layer height D dissipation parameter D(z, Z) flux distribution-function E e turbulence kinetic-energy E turbulence kinetic-energy Err MAX computational accuracy of Y-balancing F fc Coriolis coefficient F specific force F j specific force for downward motions F | specific force for upward motions Fri eddy-fraction distribution-density for downward motions Fri eddy-fraction distribution-density for upward motions F§ asymmetry specific-force due to mean-flow instabilities F T T K E specific-force due to mean-flow instabilities FjG T K E specific-force contribution due to the sub-scale shear FB buoyancy specific-force Fs shear specific-force FX" nonlocal kinematic turbulent flux of S F" 4 asymmetric specific-force ¥ s symmetric specific-force xiii 9 G gravitational acceleration H h size of turbulence domain h vertical span of V T T E H(z, Z, t) turbulent diffusivity transfer function I / turbulence invariant-measure Im dimensionless turbulence invariant-measure K k von Karman constant K eddy-diffusivity coefficient Ko eddy diffusivity at the long-time integration limit K P R O nonlocal counter-gradient flux-correction model Ki asymmetry decay parameter Ki asymmetry forcing parameter K™ T K E cascade-gain parameter KCL T K E cascade-loss parameter KF T K E forcing parameter I L scale of turbulent motions L0 stability-based length scale minimum of H Lf maximum of S L E S large-eddy simulation L F C level of free convection L N B level of neutral buoyancy L O C limit of convection M MI mixing intensity ML mixing length M M M - N C A R Mesoscale and Microscale Meteorology Division at N C A R M P mixing potential MY25 1.5 order closure with T K E prediction N rid number of layers in the turbulent domain N number of vertical layers in the finite-difference model ./V2 Brunt-Vaisala frequency N^'J^ number of daughter scales for mother [I, J] •^/(j-'^ i) number of daughters at the scale / for mother [I, J] Nsd number of layers in turbulent subdomain N C A R U.S. National Center for Atmospheric Research xiv N B L N W P neutral boundary layer numerical weather prediction P p fluctuation of mean pressure P P mean pressure P[ partial volume of downdrafts in V T T E Pr partial volume of updrafts in V T T E P B L planetary boundary layer P E potential energy PS process spectra P T M predictive transilient model R r nonlocal equivalent of RB R(z,Z) turbulent exchange function •R[n,m] discrete turbulent exchange function Rc critical Richardson-number Rc,T Richardson number including the hysteresis effect Rf flux Richardson-number Ri gradient Richardson-number RB bulk Richardson-number RB-TM parameterization of T M based on RB RPE potential-energy Richardson-number RT terminal Richardson-number Re Reynolds number Rec critical Reynolds-number s S source-layer height sgn sign function S generic source/sink terms Sw skewness of vertical velocity S B L stable boundary layer S P E shear potential-energy SPE shear potential-energy STD Spectral Turbulent Diffusivity theory S 3rd statistical moment of the vertical velocity (S-forcing) generic buoyancy forcing for c\§ (S-decay) generic decay for dS T t time te estimation time-interval for T M To turbulence time-scale Tv mean virtual temperature T3 Transilient Turbulence Theory TS" transport spectrum of S TSmr air-transport spectrum xv T K E turbulence kinetic-energy T K E - T M parameterization of T M based on diagnostic T K E T M transilient matrix T R A N version of T3 model used in comparison study T T M transilient tendency-matrix 1\S T velocity scale T turbulence kinetic-energy (T-cascadeGaingeneric turbulence energy-cascade of <9tT (T-cascadeLoss)generic decay of c\T (T-forcing) generic production of dfT by mean-flow instabilities T K E turbulence kinetic-energy U u fluctuation of zonal mean wind U u* surface friction-velocity Ui fluctuation of mean-wind component Ui U zonal mean wind Uo characteristic scaling velocity Ug zonal geostrophic wind Uh horizontal mean wind Ui i-th component of mean wind V v fluctuation of meridional mean wind V V meridional mean wind Vg meridional geostrophic wind V T T E virtual turbulent-transport eddy V T T P virtual turbulent-transport partial-volume V T T V virtual turbulent-transport velocity W w fluctuation of vertical mean wind W w nonlocal weighting factor w* Deardorff convective velocity W[ characteristic velocity for downdrafts in V T T E characteristic velocity for updrafts in V T T E W-B velocity scale w turbulent-transport velocity W vertical mean wind W2 2-nd statistical moment (variance) of vertical velocity in V T T E W3 3-rd statistical moment of vertical velocity in V T T E X Xi i-th Cartesian coordinate X nonlocal forcing function Y Y mixing potential YQ intermediate matrix in derivation of Y xvi Y r e / reference mixing potential Ypr intermediate matrix in derivation of Y Y X J J intermediate matrix in derivation of Y YLR* intermediate matrix in derivation of Y YMF intermediate matrix in derivation of Y ||Y||oo Loo norm of M P Z z height above ground ZQ surface-roughness length Zi depth of the planetary boundary layer xvii Preface The atmosphere is a fluid, so theoretically its evolution (i.e., weather) can be described by a set of mathematical equations of fluid dynamics. However, many of the physical processes in the atmosphere are not well understood. Also, to solve the equations of fluid dynamics, one must use numerical approximations. The infinite number of points in the continuous fluid of atmosphere are represented by a limited set of discrete points. Typical numerical weather-prediction (NWP) models operated by national forecast services divide the atmosphere into grid cells of air of over 20km wide and long, and a couple hundred meters deep, which are numerically represented by a single point. Anything that happens inside of this grid box is defined as a sub-grid process and is parameterized in an attempt to better approximate the true behavior of the atmosphere. Turbulence is one of those sub-grid processes that cannot be properly resolved in operational N W P models. Therefore, it must be parameterized using a turbulence model. A new class of turbulence closure - nonlocal turbulence-closure with nonlocal transport, at higher statistical order - is investigated in this dissertation. The first three chapters review the state of research relevant to this dissertation. Chapter 1 presents the properties of turbulence and turbulent transport. Typical approaches to turbulence modeling are the subject of chapter 2. Transilient Turbulence Theory (T3) is one such model, and the rationale behind choosing it as the starting point for this work is given at the end of this chapter. The details of T3 are described in Chapter 3, along with motivation for the current work. The remaining four chapters investigate whether T3 can be formulated at higher statistical order. One candidate higher-order nonlocal parameterization is given in Chapter 4. The param-eter space and properties of the new model are explored in Chapter 5. The performance of this model relative to other turbulence models is analyzed in Chapter 6. Finally, conclusions and recommendations for future work are the subject of Chapter 7. xvm Acknowledgments The author would like to express his gratitude to the following persons. To his supervisor, Dr. Roland Stull, for his patience and invaluable help that led to completion of this dissertation. To the members of his Ph.D. committee, Dr. Susan Allen, Dr. Lionel Pandolfo, and Dr. Douw Steyn, for valuable comments that contributed to the final manuscript. The author would like to thank Edward Patton and Joshua Hacker from the Mesoscale and Microscale Meteorology Division of the U.S. National Center for Atmospheric Research ( M M M -N C A R ) for providing data and software for the LES database for P B L model evaluation, used in this study for model calibration and evaluation. Last, but not least, the author would like to thank his wife, Gosia, and daughter, Marta, for years of unquestioning support and understanding. Support for this research was provided in part by the Natural Science and Engineering Re-search Council (NSERC), the Canadian Foundation for Climate and Atmospheric Science (CF-CAS) , Forestry Innovation Investment (FII), Environment Canada (EC), and the U.S. National Science Foundation (NSF). xix Chapter 1 Review of turbulence properties Turbulence is a physical phenomenon that cannot be represented exactly by numerical weather-prediction (NWP) models either today or in the near future. Two reasons are: (1) turbulence is a feature of the flow that is manifested by a certain level of 'unpredictable' disorder; and (2) the sizes of turbulent structures range from millimeters to kilometers, and are smaller than those that can be deterministically resolved in current N W P models. Hence, N W P models must rely on sub-scale parameterizations of turbulence. Turbulence is responsible for both mixing and transport of species in fluid, and is a few orders of magnitude more effective than molecular diffusion. Turbulence is second only to advection by mean wind in being mostly responsible for mixing and transporting pollutants away from the ground. Good prediction of pollutant dispersion is impossible without a reliable model of turbulent transport. To set the stage for the new turbulence parameterization that is investigated in this dis-sertation, properties of turbulence and turbulent transport are reviewed in this chapter. The specific features of turbulent transport that stimulated this work are introduced, including the mechanisms that contribute to the generation/decay of turbulence. 1.1 Proper t ies of turbulent flows 'Turbulence' is the name given to the flows that are irregular; that is, disordered in space and time. Tennekes and Lumley (1972) and Lesieur (1997) offer the following set of basic criteria to distinguish turbulence from other seemingly disorganized phenomena: • Turbulence is a property of the flow, not of the fluid; turbulent flows have certain universal characteristics, independent of the type of fluid, but depending on geometry of the flow, forcings, boundary conditions, etc. • Turbulence is a continuum phenomenon, governed by equations of fluid mechanics - but even the smallest scales of turbulent motions are much larger than molecular scales. • Turbulent flow is unpredictable in the sense of chaos theory; namely, it is sensitive to initial conditions (where small initial differences lead to large differences in evolution of the flow). The high irregularity of turbulent flow makes a deterministic approach impossible, and forces researchers to rely on statistical methods. • Turbulence is always dissipative, and its continuing existence depends on continuing produc-tion from the mean flow. Mean-flow forcings are converted into turbulence kinetic energy at 1 the larger turbulence scales, transferred through medium scales through the inertial cascade, and in the smallest scales are ultimately converted into heat by viscosity. • The mixing efficiency of turbulent flow is much larger than that of molecular diffusion. The rates of turbulent momentum, heat, and mass transfer are a few orders of magnitude higher than those of molecular transfer. • The flow must involve a wide range of spatial wavelengths. One of the most celebrated and distinct features of turbulence is its energy spectrum, filling all scales between energy injection scale to the smallest dissipation scales. • Turbulence is rotational and three-dimensional. A n eddy-like character of turbulence in-volves complex interaction of three-dimensional vortical structures (e.g.,vortex stretching) that cannot exist in two dimensional flows. (An exception is large-scale geostrophic turbu-lence governed by synoptic-scale dynamics.) 1.1.1 Pred ic t i on of turbulence The above-mentioned set of properties allows one to determine a'posteriori whether the flow was turbulent of not. However, this procedure is not sufficient to predict the onset or existence of turbulence. For prediction, one must consider the mechanisms that can generate the turbulence, such as are approximated within some key dimensionless numbers. Reynolds number Shear instability can generate turbulence. One quantity to predict the potential for shear genera-tion of turbulence against the opposing effect of viscous damping is the Reynolds number, defined as: R e = ^ (1.1) where U is a velocity scale, L is a length scale, and v is the kinematic viscosity of the fluid. The definition of U and L depends largely on the type of flow. For example, for boundary-layer flows U might be a free-stream velocity of a flow relative to the ground, and L might be the depth of the boundary layer. For non-stratified flows, experimental observations indicate that turbulent flows always occur at high Reynolds numbers. Laminar flows with common properties, such as geometry or boundary conditions, have a tendency to become turbulent at a similar value of Reynolds number, called critical Reynolds number. For example (Tennekes and Lumley, 1972) boundary-layer flows can become unstable at critical Rec ~ 600, while laminar pipe flows become turbulent at Rec ~ 2000. Thus, the Reynolds number can be used as a predictor of turbulence for specific classes of non-stratified flows. Richardson number The fact that the Reynolds number is solely related to shear instability is a major deficiency when considering atmospheric or oceanic flows. Those flows are characterized by density strat-ifications defining the effect of gravitational buoyancy force on vertical motions in the fluid. In the atmosphere, the density stratification is mostly a result of the distribution of heating/cooling due to radiative transfer, water-phase changes, and horizontal advection of air masses. In the oceans, variation of salinity is another significant factor. The effect of stratification is compli-cated because it can either supply energy to turbulence, can mitigate turbulence generation, or can prevent turbulence altogether, depending on stratification strength. 2 The static stability of the fluid can be quantified using the Brunt-Vaisala frequency, N, N2 = -9-dzQ (1.2) Q where g is the gravitational acceleration and g is the mean density of the horizontal layer of the fluid. For the atmosphere, the Brunt-Vaisala frequency is written using virtual potential temperature as a surrogate for density N2 = -§-dzev (1.3) where ft„ is the mean virtual potential temperature averaged over the horizontal plane. The Brunt-Vaisala frequency indicates { > 0 —> sub-adiabatic stratification (damps turbulence) = 0 —> adiabatic stratification (1-4) < 0 —> super-adiabatic stratification (enhances turbulence) The last stratification above is called convective instability, because the buoyancy of warm air rising and cold air sinking generates turbulence. The relative importance of shear and convective processes led to the definition of another dimensionless number: the flux Richardson number R f = wudzU + wvdzV where U and V represent horizontal components of mean wind, u, v, and w represent fluctuations of mean wind components, wOv represents mean vertical kinematic turbulent heat flux, and wu and Wu represent vertical components of mean kinematic turbulent momentum flux of U and V. At closer investigation, one can identify the denominator and numerator as the shear and buoyancy forcing terms from the tendency equation for turbulence kinetic energy (equation (3.42) introduced in Chapter 3). In a locally statically unstable atmosphere (N2 < 0 and Rf < 0), flow is unconditionally unstable as both shear and convective instability contribute to turbulence generation. In a locally statically stable atmosphere (N2 > 0 and Rf > 0), where buoyancy force opposes vertical motions, two scenarios are possible: 1) the restoring buoyancy force is weaker than shear instability (Rf < 1) and the flow will become/stay turbulent, and 2) the restoring buoyancy force is equal or stronger than shear instability (Rf > 1) and the flow will become/stay laminar. In summary ^ f < 1 —• flow tends to be turbulent (dynamically unstable) ,^ „> * \ > 1 —>• flow tends to be laminar (dynamically stable) The Richardson number may be indeterminate if its denominator is equal to zero; namely, if shear instability does not exist in the flow. In such conditions the Brunt-Vaisala frequency may be used to determine dynamic stability of the flow, which directly corresponds to static stability. In practice, the exact values of turbulent heat and momentum fluxes are not known without prior measurements or without implementing some kind of mathematical model that describes them. Therefore two approximations to the flux Richardson number are used to help predict turbulence onset. One is called the gradient Richardson number -2- r) ft R % = (dzuy + (dzvy ( L 7 ) 3 and is obtained by applying the gradient-transport model for turbulent fluxes (also known as K-theory) which relates the turbulent flux to the local mean gradient. One can also identify the numerator of i?, as the Brunt-Vaisala frequency N2; i.e., a static stability parameter. A second approximation is called the bulk Richardson number . |= Aev Az R b ~ (AU)2 + (AV)2 ( 1 ' 8 ) It is simply a finite-difference approximation of Ri for a layer of depth Az, where symbol A represents a change of a variable over distance Az, and ©„ is vertical mean of O^ in the layer. The theoretical value of the critical Richardson number for the onset of turbulence is Rc = 0.25, and was derived for inviscid heterogeneous flows by Miles (1961). In practice, the value of Rc varies depending on experimental conditions, but typically Rc ~ 0.25 is cited as an acceptable estimate. The second value, sometimes referred to as terminal Richardson number, RT, specifies conditions when turbulent flow becomes laminar and is usually specified as RT = 1. Hence the transition to turbulent flows happens when Ri < Rc, but once turbulence exists, it supports itself at higher gradient Richardson numbers up to RT, above which the buoyancy forcing prevents generation of turbulence entirely. This phenomenon is a hysteresis effect. In summary the turbulence existence criteria using the gradient Richardson number are: < Rc ~ .25 —> flow tends to be turbulent (dynamically unstable) Ri : < € (RC,RT) —> indeterminate, depends on history/triggers (1.9) > RT — 1 —• flow tends to be laminar (dynamically stable) Independent of which particular version of Richardson number one uses, none of them specifies the intensity of turbulence. Even in the extreme cases where one of the forcings dominates, such as free convection (R << —1) or forced convection (R ~ 0), only the character of turbulent forcing and some resulting properties of turbulent field, such as isotropy or homogeneity, might be deduced from Richardson number. Shear instability is always a production factor in turbulence generation, while buoyancy forcing can either intensify or suppress turbulence depending on static stability. The Richardson number is better than the Reynolds number for diagnosing turbulence in atmospheric flow. However, its dependence on local (i.e., specified at single infinitesimal horizontal layer) properties of turbulent fluxes and mean profiles limits its practical applicability to flows with simple vertical structure, such as turbulent flows with monotonic structure of mean profiles and fluxes. The commonly used gradient Richardson number also suffers from limitations of the gradient-transport approximation; namely, it is valid only for: 1) small eddies (i.e., mostly local turbulent transport) (Stull, 1991b; Arya, 1999; Pope, 2000), and 2) down-gradient local transport (Arya, 1988; Stull, 1991b; Garratt, 1992; Kaimal and Finnigan, 1994; Arya, 1999). Hence, despite including the buoyancy forcing with the shear forcing, the gradient Richardson number fails as a turbulence diagnostic tool in those real turbulent layers where positive buoyancy dominates. Namely, it fails when turbulence is driven by convective instability and characterized by existence of large eddy structures (i.e., thermals) and nonlocal turbulent transport. Interestingly, the bulk Richardson number may serve as a better indicator of turbulence existence in convective layers, provided that Az is properly specified to account for the nonlocal character of finite difference approximations. This is a feature that Stull (1984) and Zhang and Stull (1992) utilized in their nonlocal turbulent transport models, as is reviewed in Chapter 3. Nonlocal static stability The main reason behind the failure of the gradient Richardson number in diagnosis of convective turbulence is the complex structure of the mean profile of virtual potential temperature. A n ide-4 z free atmosphere + T sub-adiabatic layer 1 f adiabatic layer C B L K L super-adiabatic layer w9v Figure 1.1: A n idealized profile of mean virtual potential temperature, 9 „ , and vectors showing magnitude and direction of turbulent heat-flux, w9v, in the convective boundary layer (CBL) . z is height above ground. The turbulent C B L consists of three distinct sub-layers: a super-adiabatic layer in the lower part, an adiabatic layer in the middle, and a sub-adiabatic layer in the upper part. Positive turbulent heat-flux typically extends into the upper part of the C B L . alized virtual potential temperature profile for the convective boundary layer (CBL) is presented in Figure 1.1. Typical C B L structure consists of three distinct sub-layers: a super-adiabatic layer near the ground, an adiabatic layer in the middle, and a sub-adiabatic layer in the upper part. The observations of free-convection cases (Arya, 1988; Stull, 1991b; Garratt, 1992; Kaimal and Finnigan, 1994; Arya, 1999), show that the C B L is turbulent through its full depth with positive heat flux extending into the upper sub-adiabatic part. This convective behavior is not captured by the traditional, local, stability definitions. For example, the middle adiabatic layer is locally neutral, which according to local theory suggest that turbulence should not develop without existence of shear instability. The upper sub-adiabatic layer should be dynamically stable according to the local criteria, contradicting the observations showing both vigorous turbulent activity there as well as the counter-gradient turbulent heat flux in its lower part. Although the picture presented here is over-simplified, it nevertheless illustrates that both the local static stability criteria and the gradient-transport approximation fail to properly describe the properties of turbulence everywhere, except in the lower part of the C B L . One popular approach to overcome the deficiency of the gradient-transport approximation for turbulent heat-flux in the C B L is to add an arbitrary profile correction factor similar to the one proposed first by Deardorff (1966): ^ ~ - ( 0 z e „ - 7 C ) ( i . io) where the nonlocal correction coefficient, 7C, is estimated based on knowledge of some parameters describing surface forcing in C B L . The correction allows the gradient-transport approximation to work by creating a fictitious Qv gradient in the middle of the C B L where, in fact, there is virtually no gradient (Figure 1.1). This approach introduces a simplistic nonlocal turbulence forcing concept, according to which generation of turbulence in one part of the flow depends on the forcing in another, distant, part. 5 Stull (1991c) proposed a more general static stability criteria that works for all conditions, including convective turbulence. The approach determines static stability nonlocally, and is based on the observation that vertically displaced air parcels can move large distances due to their buoyancy. The nonlocal static stability analysis utilizes the following steps: 1. the whole profile of virtual potential temperature must be known - lack of part of the profile leads to uncertainty in stability determination, 2. conceptual air parcels are displaced adiabatically (vertically) in each direction from each point of the profile - in practice parcels need be moved only down from local minima and up from local maxima of © v , 3. continued air-parcel movement after its initial displacement should be based on its buoyancy, 4. displaced air parcels continue to move buoyantly from their origin until their paths hit a material surface or cross the profile; i.e., they reach the level of neutral buoyancy (neglecting overshoot). After tracking air-parcel movements, the stability of portions of the profile must be determined in the following order: 1. unstable: the layers where displaced parcels can continue to move under their own buoyancy, 2. stable: any sub-adiabatic layers that are not unstable, 3. neutral: any adiabatic layers that are not unstable, 4. unknown: any stable or neutral layers if parts of the profile are missing. (This category is included because observations from towers often do not span sufficient vertical range to determine static stability and turbulence potential.) The above procedure allows reliable qualitative detection of nonlocal static stability, and may be used in combination with the gradient Richardson number to identify dynamically unstable regions. Namely; all layers that are nonlocally statically unstable are also dynamically unstable, and the dynamic stability of remaining parts can be analyzed using the gradient Richardson number. For example, applying the procedure to the virtual potential temperature profile pre-sented in Figure 1.1 defines the entire C B L as being dynamically unstable, as expected. Thus one can overcome the deficiency of both the local static stability criteria and the gradient-transport approximation for convective turbulence. 1.1.2 Spectral properties of turbulence Multiplicity of scales of turbulent flow is a paramount property of turbulence. Richardson (1922) described the transfer of energy from the mean flow into turbulence at the largest scales of turbulent motions. He also noticed that turbulence energy is dissipated in the smallest scales due to molecular viscosity. The energy is continuously transferred from production to dissipation scales in a complex nonlinear energy-cascade process that fills in the whole spectrum between the energy-production and viscous-dissipation scales. The quantitative description of this Richardson cascade, for high Reynolds numbers and statistically isotropic turbulence, was first proposed by Kolmogorov (1941). He divided turbulent spectra into three distinct ranges: the energy-containing range where the energy is injected by mean flow instabilities, the inertial subrange where the energy is transferred from larger to smaller eddies, and dissipation range where energy is converted into internal energy of the fluid by 6 viscosity. Kolmogorov derived a scale-dependent form of the so-called second-order longitudinal velocity structure function for the inertial subrange, which he found to be proportional to Z 2 / 3 , where I is the scale of turbulent motions. This particular dependence implies a relationship between turbulence kinetic energy and the scale of turbulent motion in the inertial subrange, that was first derived explicitly by Obukhov (1941a,b): TKE(K) ~ e 2 / 3 K~5/3 (1.11) where K = 2ir/l is the wavenumber associated with turbulent scale Z, and e is the viscous dissipa-tion rate. The above formula was so widely confirmed by experimental results that it became a signature of well-developed turbulence. However, the span of the inertial subrange is proportional to the magnitude of the Reynolds number. The inertial subrange need not exist if the Reynolds number is not high enough. 1.1.3 Proper t ies of energy-containing eddies Even though the properties of turbulence in the inertial subrange show some universal behavior, the same cannot be said about turbulence in the energy-containing range. The structure of turbulence in this range depends strongly on the particular flow, mainly its geometry and the character of the turbulence-generation mechanism. The assumption of statistical isotropy is not valid for large eddies. Their anisotropy depends on the dominant forcing: buoyancy acts only on the vertical component of velocity perturbations, while shear terms affect mainly horizontal components of velocity perturbations. For positive Richardson numbers, the eddies are compressed in the vertical due to both the damping by buoyancy and the anisotropy of shear forcing. Flows with large negative Richardson number are characterized by vertically elongated structures forced by convective instability, with the amount of anisotropy corresponding to the magnitude of Richardson number. The anisotropy tends to diminish as the magnitude of the Richardson number decreases, although due to anisotropy of shear forcing one cannot expect neutral layers to be fully isotropic. The energy-containing eddies often exhibit evidence of well organized, i.e. coherent, struc-tures. Coherent structures are most visible for free convection, where well-organized positively-buoyant thermals develop throughout most of turbulent layer, and are accompanied by compen-sating downdrafts. The strong coherency of thermals allows them to overshoot and penetrate deeply into the overlying stable layer, that caps the mixed layer, and forces air from that layer to be mixed down into the convectively unstable part, forming a relatively deep entrainment zone at the interface between the unstable and stable layers. Large horizontally-dominant eddies can also develop in stable layers, although they do not show the same level of coherence as convective eddies. The scales of energy-containing eddies also depend on the stability of the turbulent layer. Kaimal et al. (1976) analyzed turbulence spectra in the convective boundary layer, and found that the wavelengths of peak spectral energy for horizontal wind components are independent of height, and are equal to 1.5ZJ (where Zi is the depth of boundary layer defined as the height of the capping temperature inversion). The peak spectral amplitudes also grow slowly with increasing height. For vertical velocity, the wavelength of the spectral peak grows exponentially toward a value of 1.5zi with height until 0.5ZJ, and stays approximately constant through the upper half of the convective boundary layer. The peak spectral amplitude grows rapidly from the ground up to 0. 1ZJ , and increases relatively slowly in the rest of the boundary layer. For temperature spectra, the wavelength of the spectral peak shows similar behavior, but approaches 1.5ZJ at a lower level; 1. e., at O.IZJ . The peak spectral amplitude however behaves differently, as it drops steadily until 7 0.5zi, then stays constant up to 0.7Zi where it rises again due to entrainment processes in the upper convective boundary layer. This behavior of the vertical velocity and temperature spectral peaks on z in the lower part of the convective boundary layer reflects the influence of surface boundary conditions. Namely, vertical velocity must become zero at the solid surface. Nevertheless, those results indicate that the length-scales of energy-containing eddies in the C B L depend primarily on the coherent convective structures extending throughout most of turbulent layer, and hence they scale well with z%. The nonlocal character of convective buoyancy forcing and strong coherency of large-scale eddies in convective layers lead to another feature of large eddies; namely, a strong asymmetry of vertical motions. The asymmetry can be represented by nonzero skewness of vertical velocity S« = =375 d-12) which is typically close to 0.6 (Wyngaard and Weil, 1990) in C B L s . In C B L s , positive skew-ness implies that updrafts have higher vertical velocity and cover smaller horizontal areas than downdrafts (Greenhut and Khalsa, 1982). In convective layers decoupled from the surface, the skewness may depend on the relative location of the dominant forcing; e.g., a cloud layer cooled at the top due to radiative cooling might have a negative skewness associated with faster downdrafts covering smaller horizontal area (Cotton and Anthes, 1989). 1.2 Proper t ies of turbulent t ransport 1.2.1 General properties The dispersion of any quantity in a turbulent flow is directly related to the properties of tur-bulence. As emphasized by Stull (1993), turbulent transport, and hence turbulent dispersion, originates from the advection terms in the governing equations. The advective origins of turbu-lent dispersion imply that its properties are strongly affected by advective features of turbulent motions including: 1. turbulent transport of any conservative quantity can be related to the transport of air mass, 2. turbulent transport is realized by the cumulative action of eddies of different scales and of different levels of coherency, 3. the vortical character of turbulence, and mass conservation, imply that parcels moving from one level to another will be accompanied by a return flow in the opposite direction. 4. turbulent transport in the atmosphere is often statistically anisotropic and inhomogeneous. Those features clearly distinguish turbulent dispersion from molecular diffusion. Despite that, turbulent dispersion is often approximated as a diffusion-like process (K-theory). The dispersive properties at different stabilities are best illustrated by following the vertical spread of a passive tracer, released from a point source such as a smoke-stack (Arya, 1999). In stable conditions, the plume away from the source has a near-Gaussian vertical distribution (fanning). The vertical extent of the plume increases away from the source at a rate determined by the strength of stability; a stronger stability implies a narrower vertical distribution of pollutant, as it has stronger damping effect on energy-containing eddies. The center of the mass of this passive tracer stays at a constant level, provided that mean motions on the scales larger thaii 8 turbulence do not affect its location. In unstable conditions, the coherency of large eddies leads to looping; i.e., a situation where the emissions are injected into traveling updrafts and downdrafts, and hence the level of the plume local center of mass varies accordingly. Far away from the source, a well-mixed layer of pollutant is created filling the full depth of the C B L . The vertical stability structure of the flow can influence not only the effectiveness of turbu-lent transport but also its homogeneity. First, turbulent transport is limited to the depth of the turbulent layer. Strong stable zones outside of this layer can entirely suppress turbulence generation and turbulent dispersion. This so-called trapping effect is most dominant in C B L s topped by a strong capping inversion. Vertical variation in the strength of stability can also lead to strong inhomogeneity of turbulent transport in stable layers. Lofting is an extreme case where a neutral layer lies on the top of a strongly stable layer, and the tracer is released in-between those layers. The tracer spreads nearly exclusively into the neutral layer. Overall, the variations in vertical structure of stability imply anisotropy and inhomogeneity of turbulent transport, and hence non-Gaussian vertical dispersion. The spectral properties of turbulence also influence the structure of turbulent transport. The inertial-subrange scaling of turbulence implies that the large scales have stronger impact on turbulent transport than small isotropic scales. As turbulent transport is determined mostly by those large scales, one can expect that their specific properties are to be reflected in the overall turbulent dispersion. In stable boundary layers (Caughey et al., 1979), turbulent fluxes are strongest near the surface (in the layer of strongest shear, in agreement with highest amplitudes of turbulence spectra) and decrease steadily towards the top of boundary layer as do the amplitudes of the turbulence spectra. In CBLs , Greenhut and Khalsa (1982) determined that throughout most of the layer the well-defined coherent updrafts and downdrafts are responsible for nearly 90% of turbulent fluxes, even though they occupy only 40% of the horizontal area. 1.2.2 Specific features of convective transport Mahrt (1976) observed that the potential temperature is usually well mixed in the middle part of C B L s , while the mean specific humidity is often not. He suggested that, among other factors, an entrainment process at the top of the C B L might be responsible for that. Wyngaard and Brost (1984) confirmed that hypothesis, and further attributed that behavior to vertical asymmetry of buoyancy forcing. They introduced the idea of separating the C B L diffusion processes into two distinct components: 'bottom-up' diffusion driven by the surface flux, and 'top-down' diffusion driven by entrainment-induced flux at the top. Wyngaard and Brost (1984) and Moeng and Wyn-gaard (1984) proposed formulas for eddy diffusivity that scaled with z,, and which described the strong asymmetry between top-down and bottom-up processes and the counter-gradient bottom-up transport in the upper part of the C B L . Numerous later studies (Wyngaard, 1987; Sawford and Guest, 1987; Weil, 1990; Wyngaard and Weil, 1990; Holtslag and Moeng, 1991; Cuijpers and Holtslag, 1998) have shown that three features of C B L turbulence influence transport asymmetry: 1) entrainment at the top of C B L due to penetrative convection, 2) skewness of vertical veloc-ity fluctuations (Equation (1.12)), and 3) vertical inhomogeneity of variance of vertical velocity fluctuations w2. Coherency and skewness of large eddies in C B L s cause another peculiar turbulent transport property; namely, the level of maximum concentration of the tracer plume away from the source depends on the height of the tracer source. Willis and Deardorff (1976, 1978, 1981) observed this non-Gaussian behavior in laboratory tank experiments. For near-surface releases, the maximum-concentration level stays initially near the ground, and then lifts off rapidly to 0.8ZJ and gradually descends to 0.5ZJ. In the cases when passive tracer is released from elevated sources in the bottom half of the C B L , the mean concentration level descends steadily to the surface, then rises 9 above 0.5ZJ, and finally descends back to 0.5zi. Similar features were observed in the large-eddy simulations of Ebert et al. (1989), and in the real atmosphere (Eberhard et al., 1988). The above features of convective turbulent transport result from specific characteristics of convective forcing, and can be attributed to the asymmetry of buoyancy forcing that is not captured by simple gradient-transport models. A more complex turbulence model should account for the features of large-scale turbulent eddies, such as skewness and coherency, as well as vertical inhomogeneity of turbulence. This is one of the motivations for this dissertation research. 10 Chapter 2 Review of turbulence modeling 2.1 T h e turbulence closure problem The generic prognostic equations for dry incompressible turbulent flow (Stull, 1991b) are (in Einstein's summation notation): where 0 is mean potential temperature with fluctuating part 6, Ui is mean wind component with fluctuating part u^, C is a passive tracer with fluctuating part c, and S are generic mean source terms for those variables. A l l of the variables depend on time and location. On the left-hand-side are the tendency and advection terms. On the right-hand-side are the net source terms (e.g., radiative heating, buoyancy, Coriolis force, etc.) and the turbulent transport (turbulent advection) terms. The last terms contain covariances (i.e., second-order statistics), which are unknown. Forecast equations for each of these second-order terms can be derived, but they contain nine or more new third-order terms. Repeating this process ad infinitum, one finds that either an infinite number of equations is needed to describe turbulent flows, or if one stops with a finite number of equations there will always be some unknown statistical terms. This difficulty is known as the closure problem of turbulence, and has been recognized during the past 100 years as one of major unsolved problems of classical physics. 2.2 R e v i e w of turbulence closures 2.2.1 Overv iew One practical way to work around the closure problem is to retain only a small number of forecast equations, and parameterize the remaining unknowns. For example, one can retain the forecast equations for only mean variables, and parameterize the unknown second-order covariances as functions of known mean quantities and of independent variables such as height. The result is called first-order closure. While third-order closure is the highest statistical closure that has appeared in the literature, the exponentially-increasing complexity yields diminishing returns in improvements of forecast accuracy of the lower-order terms such as mean temperature. For this reason, first-order closure is most widely used. ^ 0 + 1 ^ , 0 dtUi + UjdXjUi dtC + UjdcUi SQ - dXjUj6 Su - dXjUjUi Sc - dXjUjC (2.1) (2.2) (2.3) 11 Table 2.1: Summary of turbulence closures (after Stull (1994b)). 6 denotes mean potential temperature, 9 denotes perturbation of mean potential temperature, w denotes perturbation of mean vertical velocity, dt denotes time derivative, and dz denotes height derivative. Economy Detail Nonlocalness (N) wi N2 none or local mixing bulk or nonlocal cor-rection 2-stream full nonlocal Statistical order (S) symbol sample equations 1 i i So G = param. or dzS = param. similarity theory (e.g. log wind profile ) bulk or slab mass-flux Si dtB = . . . - dxw8 •wB = param. K-theory modified gradient top-down/ bottom-up spectral diffusivity transilient integral turbulent adjust-ment mixing-length the-ory horizontal roll 2-stream direct interaction approx. orthonormal ex-pansion Si.5 dtB = . . . - dzw6 dt(TKE) = param. w6 = param. T K E k-c Yamada level 2.5 higher-order nonlocal (this work) s2 dtQ = . . . - dzw6 dt^9 - . .. - dzw2e w^B — param. 2nd-order closure 2nd-order mass flux Sa dte = . . . - dzwe atwd = . . . - ozw2e dtw2e - . . . - dzw'^e w39 = param. 3rd-order closure Within any statistical class of closure, a wide variety of parameterizations are possible for the unknown covariances. The covariance of velocity with any other variable can be interpreted physically as a turbulent flux of that variable. For example using an overbar to denote a mean value, and lower-case letters to denote perturbation (i.e., a deviation from the mean) the covari-ance of vertical velocity perturbation (w) and potential temperature perturbation (6) is w9, and represents the vertical heat flux in kinematic form. One form of first-order closure is to assume the heat flux flows down the local gradient of the mean potential temperature (0): w6 = —KdzQ. This approach is called K-theory, named after the symbol K used to represent the eddy thermal diffusivity. It is a local closure because the heat flux at any height (z) depends on only the local temperature gradient at that same height. Nonlocal closures are also possible, where the heat flux at any height depends on the potential temperatures at many heights. In the atmosphere, large coherent eddy structures such as thermals of warm rising air can transport heat nonlocally across large vertical distances, in a process that is poorly described by K-theory. For this reason, nonlocal closures have gained acceptance in the past decade. One such nonlocal closure is called Transilient Turbulence Theory (T3). Table 2.1 from Stull (1994b) shows the various orders of statistical and nonlocal closures that have appeared in the literature. Most turbulence parameterizations consist of two parts, a framework, and a description of the free parameters. For K-theory, the framework is the local down-gradient equation w9(z,t) = —K(z,t)dz9(z,t). The free parameter K can be modeled as a function of height, Richardson number, wind speed, turbulence kinetic energy (TKE) , etc. For T3 in discrete (finite difference) form, the framework is a matrix representation of nonlocal mixing 0[j](£ + At) — 2~2jLi ^[ij]^) w n e r e i n * n ^ s discrete formulation, the index i or j 12 represents the height of the center of the layer in a column of grid points (e.g., = iAz — Az/2 for equal-height layers), [©] vector represents potential temperature (but could be any variable that is conserved under adiabatic vertical movement), and [C] is a matrix of mixing coefficients that is known as a transilient matrix (TM). Each element C[jj] of the T M indicates what fraction of air ending at destination cell i came from source cell j. The [C] matrix is also known as a transfer matrix or source-receptor matrix. The parameters can be modeled as a function of height, Richardson number, wind speed, T K E , etc. Since T3 was devised about 20 years ago (see review by Stull (1993)), most of work has been on first-order statistical closure. Namely, the C^j j parameters were modeled as a function of vertical differences of mean variables (wind, potential temperature, etc.) at height indices j and i. Also, for simplicity, the [C] matrix was assumed to be symmetric; namely, the amount of air mixing from j to i was assumed to equal the amount from i to j. However, measurements of the T3 matrix for convective boundary layers in a lab tank (Cenedese and Querzoli, 1997) and numerical large-eddy simulations (Ebert et al., 1989) suggest that the matrix becomes asymmetric for free convection. Also, the change of this matrix for different-duration time-steps cannot be fully explained by the simple assumption that the matrix for 2At just equals the matrix for lAt applied twice (Ebert et al., 1989). These parameterization deficiencies might be alleviated by utilizing a higher-order statistical order for this nonlocal closure. This is the overall goal of this dissertation. 2.2.2 Local closures The most common local closure is called the gradient-transport model, or simply K-theory. It was developed by analogy to molecular diffusivity. It is a first-order statistical closure, where the turbulent flux of any mean quantity E, with perturbation £, is expressed as a function of its mean gradient, and of a parameter called the eddy-diffusivity coefficient K: x uJ^-K&Z (2.4) In the simplest formulation K is constant. A more practical extension of K-theory is the mixing-length model. It goes a step further and parameterizes AT as a function of absolute shear \djUi\ for i / j, by introducing a mixing length parameter, I. For example, in the atmospheric boundary layer one can parameterize K as K = l2\dzUh\ (2.5) where Uh is a horizontal velocity component. Mixing length may be expressed in different ways. The simplest one is I2 = k2z2 (2.6) where k ~ 0.4 is the von Karman constant. K-theory and its extensions are often described as small-eddy closures, because they fail when large eddies (such as thermals) are present in the real flow. Although some more sophisticated parameterizations of K and I attempt to account for the existence of convective instability, there is considerable evidence (Stull, 1993) that they fail in various situations, especially where counter-gradient turbulent transport is observed. One-and-a-half order closure (e.g., Mellor and Yamada (1974, 1977)) is presented here as an example of higher-order closure. The key difference between this closure and K-theory is in carrying prognostic equations for variance of velocity fluctuations (turbulence kinetic energy 13 e = \uf) and potential temperature (92). The governing equations are defined as follows: = -K,:(e,62)dzZ (2.7) dte = —wudzll — wvdzV + — dzw(p/g + e) — e (2.8) = -2w9dzS - dzw92 - 2se (2.9) = \Kel'2dze (2.10) we2 = Keel'2dz¥ (2.11) e3/2 e (2.12) (2.13) K where E represents any prognostic mean quantity with perturbation £, e are dissipation rates for the respective quantity, and A are empirically obtained length-scale parameters. One may see that fluxes of e and 92 are still expressed in the K-theory fashion; i.e., by the gradients of those quantities and a proportionality coefficient. This is a general trend in the local higher-order modeling of turbulence. Summarizing, based on this simple example one can see certain features of local closures. Higher-order closures will demand more computations than lower-order. Even though closure as simple as one-and-a-half usually improves the prediction of shear generated and convective structures, it still suffers from the inability to predict counter-gradient and zero-gradient transport of any quantity, because the direction of the local flux is still modeled based on the down-gradient transport assumption. (Notice that K would need to take unphysical values for some flows, such as negative for counter-gradient transport, or even infinity for profiles with zero gradients.) Hence, turbulent transport due to large eddies is still poorly simulated. Finally, those schemes are not able to simulate nonlocal turbulent advection in any sense, such as by convective thermals or other large coherent structures. This poses severe constraints on integration time-step of any numerical model, since any quantity may be moved only to the nearest grid during a single step. 2.2.3 Nonlocal closures Nonlocal turbulence closure models (Stull, 1993) represent the combined action of all-size eddies, from smaller turbulent eddies to larger coherent convective structures. These large eddies are generated by nonlocal properties of the flow, and result in nonlocal turbulent transport. One may classify nonlocal closures in three basic ways (Stull, 1993). First, by the statistical order of the closure, as in local closures. Second, by the form: analytical or finite-difference'; where analytical closures assume different types of integro-differential forms, and finite-difference forms are represented by discrete matrix operators. Third, by the underlying assumptions about the nature and structure of turbulence: a priori vs. responsive. The a priori models are based on a prescribed turbulence spectrum obtained from numerical and field experiments. The respon-sive models relate turbulent transport to actual properties of the flow, such as the tendencies, gradients, and turbulence intensity. A priori models work fairly well in idealized cases, but are rarely applicable to more realistic non-stationary and non-homogeneous flows. Overview of Transilient Turbulence Theory ( T 3 ) A complete description of T3 is given in Chapter 3. Here, only the basic elements of the finite-difference model are presented as an example, to allow comparison with local models and with 14 other nonlocal turbulence closures (see Appendix A) . The basic framework of T3 (Stull, 1984, 1993), as was presented earlier, is that the evolution of any conserved quantity, E , during a finite time interval, A t , is N H [ t ] (t + A t ) = £ c M ( t , At)Sfl(t) (2.14) where N is the number of layers in a discrete model. Recall that [C](t, At ) is called the transilient matrix (TM). In existing parameterizations, computation of the T M consists of two major steps: (1) defining the mixing potential (MP), and (2) normalizing the M P into a T M to ensure mass conservation. The derivation (Stull and Driedonks, 1987) and the recent parameterization (Inclan et al., 1996) of elements of the M P matrix are based on a nonlocal diagnostic form of the tendency equa-tion for turbulence kinetic energy ( T K E ) . This particular equation is used because it describes explicitly the state of turbulence due to two mean-flow instabilities: shear instability and con-vective instability. In essence, off-diagonal elements of the M P matrix, [Y](t), are proportional to corresponding elements of the nonlocal T K E Y M ( I ) ~ (irk) T K E M ( * ) I F * ^ J (2-15) [TKE] (t) is a function of nonlocal mean-flow instabilities defined in terms of nonlocal differences of mean-flow properties T K E M ( t ) = / (A4 / ( t ) , A ^ V ( t ) , AMe„( t ) , A^z, Rc, g, D) (2.16) where Rc is critical Richardson number, g is Earth's gravitational acceleration, D is dissipation parameter, and A l J ' ( ) represents nonlocal difference, e.g. A}'h = z[i]—z[j]. The physical constraint for kinetic energy, which requires it to be always positive, sets the condition TKE[y] = max(0,TIKE [ i ) i ] ) (2.17) and limits the range of turbulent transport to the layers of the models where diagnosed T K E is larger than zero. A n additional renormalization procedure is applied to the M P matrix to mimic the theoretical properties of the T M . Namely, the T M should approach the identity matrix for infinitesimal time-steps, and it should approach the well-mixed matrix for infinite time-steps. Therefore, the off-diagonal elements of [Y] are recalculated according to the formula: Yh= J^V1 if 3 (2-18) where 8 is a parameter determined from numerical experiments. Finally the diagonal elements of the M P matrix, which represent subgrid-scale turbulence, are defined in the following way: Y M = Y r e / (2.19) where Y r e / is also found from numerical experiments. Note that due to the definition of [TKE] for this particular T3 parameterization, the M P matrix is symmetric with respect to its indices; i.e., Y M = Y[jA. The procedure sketched above yields the M P matrix. However, state and mass conservation requires that transilient matrices are doubly-stochastic; i.e., each row and column of the matrix 15 must sum to 1. Hence, the last step of building a T M from the M P involves ensuring mass conservation properties of T M . Two alternative algorithms were devised to achieve the goal, as will be presented in Chapter 3. In summary, T3 is not a turbulent diffusion concept, but describes turbulent dispersion as an advective-like process. It parameterizes nonlocal turbulent advection in physical space as a function of nonlocal mean properties of the flow in physical space; i.e., it is a responsive parameterization. It is not limited to homogeneous turbulence, and it simulates the combined action of all eddy sizes from the smaller chaotic turbulent eddies to the larger coherent structures. Nevertheless, if no large eddies exist in the flow, then T3 mimics the small-eddy transport of K -theory. (For example, one can expect that T3 will approximate the logarithmic wind profile if the surface layer is sufficiently resolved, although it is not a typical application of a turbulence-closure scheme in atmospheric models.) T3 is also a first-order statistical closure, and it can be expressed in either analytical or algebraic form. Finally, it is nonlocal in time but still allows even the weakest eddies to contribute to the turbulent transport across the entire turbulence domain. Advantages of T3 over other nonlocal closures There are a few factors that make T3 unique, compared with the other nonlocal closures described in Appendix A , even though it is still a first-order closure. First of all, T3 is a responsive parameterization, taking into account the real strength of both common sources of flow instabilities; i.e., shear and/or convective instabilities. Also, it does not rely exclusively on turbulence scales, or on the turbulence energy spectrum. Hence, it does not prescribe any universal spectrum of turbulent eddies like other nonlocal models. In addition, unlike the other nonlocal closures, T3 does not rely on any external information, such as surface similarity scaling, about the turbulence characteristics. So, it is not only a responsive but also a self-contained parameterization. T3 explicitly simulates nonlocal turbulent fluxes in the form of nonlocal turbulent advection rather than diffusion. T3 is the only nonlocal closure that explicitly models the structure of turbulence, even though it does not follow the evolution of the turbulence. It explicitly simulates all the resolved structures; i.e., transilient matrix is a function of both location, and the distance between non-neighboring points, and may be applied to homogeneous as well as heterogeneous turbulent flow, and to finite-range turbulent exchange. T3 is the only nonlocal closure that attempts to modify the nonlocal flux structure in order to account for varying integration time-step. This limits the problem of unrealistically large rate of nonlocal turbulent motions for infinitesimal time-steps. T3 in discrete form is nonlocal in time. Time-step dependence is carried in the parameteri-zation of the transilient matrix. This makes it more flexible in numerical applications, although time-step dependence of the transilient matrix is a complex, and not yet entirely solved, prob-lem. Because T3 models explicitly nonlocal turbulent advection, it is more suitable for numerical modeling of turbulent-scale convective problems, rather than deep thunderstorm convection. A l l of these features make T3 a very interesting topic, for both applications and further research and development. It is recognized in the scientific community as an alternative to local closures in numerical modeling. As a result, there is a growing demand for further research both in the area of applications, and in theoretical development of the T M parameterization. The next chapter will expand the review of T3 beyond the brief overview that was presented so far. 16 Chapter 3 Review of Transilient Turbulence The foundation of Transilient Turbulence Theory (T3) and a review of parameterizations for the transilient matrix (TM) are presented in more detail in this chapter. There is a significant amount of literature devoted to T3 and its applications. The following short review is by no means complete, but it summarizes the most important works in the areas of development of the theory and its applications. Stull (1984) introduced the concept of eddy-mixing across finite distances and defined the fundamental properties of T3. The article also presented the first parameterization of the T M as a function of Richardson number. Stull and Hagesava (1984) compared properties of this param-eterization to a turbulent adjustment scheme. The dispersion properties of an analytical form of the T M were discussed by Stull (1986). Those three works form an introductory foundation for T3. The further development of the T3 parameterization can be traced throughout the follow-ing works: Stull and Driedonks (1987) derived another parameterization for the T M based on a diagnostic equation for the T K E , and they compared the resulting simulations with field obser-vations. Stull (1987) described a full numerical algorithm for a diagnostic-TKE parameterization of the T M , including F O R T R A N source-code, ready to implement in numerical models. Stull (1991b, 1993) presented an updated version of T3, including modifications necessary for non-equally spaced vertical grids, and compared it with other nonlocal closures. Zhang and Stull (1992) proposed another parameterization for the T M based on the bulk Richardson number, and compared its performance with a diagnostic-TKE version of the T M . Inclan et al. (1996) presented the most recent version of diagnostic-TKE version to parameterize the T M . The first attempts to formulate a symmetric parameterization of T M based on a prognostic equation for T K E were presented by Modzelewski and Stull (1998a,b). Switching from parameterization to application, the following works exemplify utilization of T3 in different areas of numerical modeling. Stull and Kraus (1987) applied T3 to the simulation of turbulent transport in the upper ocean. Raymond and Stull (1990) applied T3 in a mesoscale N W P model. Inclan et al. (1996) applied T3 to study turbulent transport in forest canopies and, Modzelewski and Stull (1997) applied T3 in an air-pollution dispersion model for passive scalars. Multiple numerical experiments were performed to analyze the behavior of T3 in different forc-ing conditions. Stull (1988) discussed turbulent dispersion associated with idealized symmetric and asymmetric transilient matrices. Ebert et al. (1989) analyzed transilient matrices calculated from L E S experiments. Stull (1991a) compared LES-estimated vs. parameterized transilient ma-trices. Stull and Bartnicki (2000) used an air-pollution dispersion model to analyze dispersion properties of T3 based on diagnostic T K E . Although Stull (1984) developed an extensive formal procedure leading to a continuous repre-sentation (using integro-differential equations) of generalized nonlocal flux, it is intended here to 17 concentrate mainly on the discrete form of the parameterization. Hence, throughout this chapter, the finite-difference view of a fluid represented by N layers is assumed, with a discretized system of equations in both space and time. The following subset of variables is used to describe the state of the mean flow: wind components [U] and [V], virtual temperature [Tv], and virtual potential temperature [0„]. The variables are defined as area averages of horizontal layers of thickness [Az], each layer centered at [z]. (For the special case of uniformly spaced grid, Az without a subscript is used to denote layer thickness.) For example, 0„[j] represents the mean virtual po-tential temperature of the layer i, which has thickness Az^y and which is centered at height z^y Nonlocal differences are expressed by means of a nonlocal difference operator A l J ; e.g., distance (z[jj — z^j) is expressed as Al,Jz. Nonlocal averages are expressed by means of a nonlocal average operator A t J ; e.g., nonlocal average of virtual potential temperature from two layers, ^e"[ ,)^e"|j)^ is expressed as A l ' J 0„. A i is used to denote the integration time-step of the numerical model. 3.1 Fundamentals of Transil ient Turbulence Theo ry The recognition of four basic facts about the nature of turbulence forms the foundation of T3. The general properties of turbulence have already been discussed in Chapter 1. Hence, only the features of turbulence that have a direct implication on formulating T3 are reiterated here: 1. Turbulent transport has an advective character. Namely, the turbulent-flux terms in the tendency equations for different quantities originate from advection terms in the correspond-ing initial governing equations. Hence, the traditional treatment of turbulent transport as a diffusive process, although possible in limited circumstances, may alter the modeled trans-port properties in inappropriate, unphysical way. 2. Turbulence has a multi-scale character. Turbulent eddies of different sizes, spanning from small-scale, energy-cascade-driven vortices, to large-scale, convectively-forced coherent struc-tures, contribute to cumulative effect of turbulent flux. Hence, turbulent transport has a multi-scale character as well. 3. The generation of turbulence depends on nonlocal properties of the flow. The classical view, local in character, does not take into account the fact that turbulent motions can be initiated elsewhere in the flow and can span across sub-adiabatic or shearless regions of the flow. Hence, turbulence generation is affected by nonlocal forcing. 4. The turbulent transport of all quantities is caused by the turbulent advection of fluid mass. Hence, if one can simulate the vertical turbulent transport of mass in the fluid, than one can use the same expression to simulate the transport of any quantity that is conserved during the vertical motion of fluid. Namely, when two different air parcels mix, the relative masses of these two parcels control the final mixture of all properties including moisture, heat, pollutants, and momentum. Other than the properties of turbulence itself, one must also take into account the numerical constraints of the modeling framework, within which the turbulence closure is to be applied. Figure 3.1 represents a vertical cross-section of single column in the grid of a typical medium-resolution mesoscale weather-forecast-model with horizontal sizes of grid cells, A x ~ 20km, and typical vertical span of grid cells, Az ~ 50m, near the bottom of the model. The typical time constraints of model integration are based mainly on mean advective winds and the integration time-step for such a setup, which is typically over 100s. The thick-line boxes represent approx-imately the horizontal and vertical span of isotropic circulations with horizontal and vertical 18 •Ax ~ 20km- w ~ 4 ^ 9 8 ™ ~ 3 f 7 r ™ ~ 2 f 6 5 4 3 2 1 L A z ~ 50m A t ~ 100s Figure 3.1: Vertical cross-section of a single column of a typical medium-resolution mesoscale-forecast-model grid. The thick-line boxes represent scopes of circulation of typical convective eddies (with horizontal to vertical aspect ratio of roughly 2) with various velocity scales as indi-cated by w, preserving approximately the aspect ratio of grid sizes. velocities ranging from 1 to 4™- One can see that even the horizontal span of the eddies is easily contained within the column, except on the very boundaries. However, the circulations penetrate through multiple layers in the vertical direction. Hence, due to the anisotropy of the horizontal vs. vertical resolutions in the models, turbulent eddies are not resolved horizontally, even though they span over multiple cells and are resolved in the vertical. The modeling framework not only sets the requirements on the structure of turbulence closure, but it also limits its applicability. It is not appropriate to use a turbulence closure that duplicates the effects of resolved scales. If the horizontal size of the grid cell in Figure 3.1 was reduced to equal the size of its vertical thickness, and if the integration time-step reduced accordingly, then the circulations would be directly resolved in both directions, and nonlocal turbulent-transport approach would not be needed. In summary, a realistic turbulence closure should take into account properties of turbulence, properties of the turbulence forcing, and properties of the modeling framework within which it is utilized. 3.1.1 Transilient matrix concept Taking the above considerations into account, Stull (1984) formulated a nonlocal turbulence closure model applicable to numerical models with anisotropy of grid resolution in the horizontal versus the vertical direction, and named it Transilient Turbulence Theory, or T3 in short. The word 'transilient' originates from the latin verb 'transilire' meaning 'jump over' or 'leap across'. The structure of T3 is based on the fact that the turbulent transport can occur not only between neighboring layers but also 'directly' between those layers separated in space. This nonlocal transport takes place due to both numerical constraints of the modeling framework and the multi-scale character of turbulent transport. The framework for T3 was already given in Chapter 2 (see Equation (2.14)), and is rewritten 19 C[4,5] r C[i,5] C[l,4] layer 5 C[4,4] J J layer 4 layer 3 layer 2 layer 1 Figure 3.2: Graphical representation of the transilient coefficients for a 5-layer model. The coef-ficients define the relative contribution of eddies of different sizes to the total turbulent transport between different layers of the model. here to promote further discussion: %](* +At) = £ e M ( i , Ai )S y ] ( i ) (3.1) The physical limitation of the Equation (3.1) is thoroughly discussed by Larson (1999). He links the T M to the Green's function of the advection-diffusion equation, and points out that the information about the sub-grid distribution of a transported species is lost every time that Equa-tion (3.1) is applied. This particular feature, commonly referred to as the convective-structure memory problem, is commonly shared by all closures based on the mean properties of the flow, and cannot be solved without rising the tensor-order of Equation (3.1) to a higher level. Such an extension, although theoretically possible within the T3 framework, would significantly increase its computational cost and could make it unfeasible for many practical modeling applications. Each element of the T M , C[ij], called a transilient coefficient, is a dimensionless number representing the portion of the mass of the destination layer i (first index) that came from the source layer j (second index) during time A t . This destination-oriented definition was introduced by Raymond and Stull (1990), as a modification to original definition (Stull, 1984), in order to allow for non-uniform masses of layers of numerical model, resulting from their variable height as well as variable density of the fluid. The meaning of the transilient coefficients is illustrated in Figure 3.2. The first index of the transilient coefficient always represents the destination of turbulent mixing, and the second index represents its source. For example, coefficient <C[5^ represents the fraction of mass in layer 5 that was transported from source layer 1, and coefficient represents the fraction of mass in layer 1 that was transported from layer 4. The elements with the destination index larger than the source index correspond to upward motions (sketched in Figure 3.2 on the right side of the grid cells), while the elements with the destination index smaller than the source index correspond to downward motions (illustrated on the left side of the grid cells). The diagonal elements of the T M , such as C ^ j , correspond to unresolved sub-grid-scale turbulent transport (i.e., unresolved in both the vertical and horizontal directions; namely, air remaining in the same 20 D5 0.8 0.2 0 0 0 D4 0.2 0.8 0 0 0 D3 ( o 0 1 0 o ) D2 0 0 0 0.6 0.4 D I 0 0 0 0.4 0.6 S5 S4 S3 S2 SI \ Figure 3.3: Example of a transilient matrix for a 54ayer model with two turbulent sub-domains separated by a non-turbulent layer. ' D ' indicates destination index and 'S' indicates source index. A n arrow in the lower-right corner indicates the main diagonal of the matrix. The turbulent subdomains are indicated by thick-line boxes. A non-turbulent domain is represented by an oval-finished cross. The matrix is flipped, as compared to the typical mathematical notation, so that the main diagonal has the same direction, but the lowest indices are in lower-right corner. The notation was introduced by Stull (2000) for easier interpretation of a Transilient Matrix; namely, higher in the matrix corresponds to higher in the atmosphere. grid cell). In summary: { i > j => upward motions i < j downward motions (3.2) i = j unresolved sub-grid-scale motions The difference between the indices of any transilient coefficient corresponds to the turbulent-transport range. The further the element is from the diagonal, the larger transport jange it represents. For example, the absolute value of the difference between indices of the transilient coefficients is proportional to transport range for uniform vertical resolution: \i — j\ oc transport range (3.3) Thus, coefficient C[ 4 2 ] represents shorter range of nonlocal turbulent transport than coefficient C[ 1 5 ] . Coefficient C[ 4 i 5] represents local transport between two neighboring layers 5 and 4. The T M can also indicate which subdomains are not turbulent (i.e., laminar), and can indicate the intensity of turbulence elsewhere. For any fixed A i , the values of the off-diagonal transilient coefficients represent the relative intensity of the mixing process. The larger an off-diagonal transilient coefficient, compared to the others, the more intense is the mixing it represents. Non-turbulent layers are represented by diagonal elements of the matrix equal to 1 and all other elements of the corresponding rows and columns equal to zero. In summary: ^ f = l =» non-turbulent layer , . VA ' | < i =^> p a r t of a turbulent subdomain In the case when fluid comprises a mixture of still and turbulent layers, each turbulent sub-domain is represented by a sub-matrix of full T M . For example, the T M in Figure 3.3 represents 21 two turbulent sub-domains separated by a non-turbulent layer. For this illustration, the turbulent mixing in the sub-domain of layers 1-2 is more intense than the mixing in the sub-domain of layers 4-5. Physical constraints of the transilient matrix (TM) Since the T M represents turbulent transport during a finite time interval, its elements should reflect the physical properties of turbulent mixing during that interval. However, the transilient coefficients and resulting T M must also satisfy a set of physical constraints. First, transilient coefficients represent the fraction of air ending in the destination layer that came from a source layer. Turbulence is always associated with increase of entropy and random-ness, therefore the process of un-mixing is neither possible nor physically realistic. Therefore transilient coefficients cannot be negative. Also, for quasi-incompressible flow, it is not possible to have more air coming into a grid box than it can hold, implying that any transilient coefficient cannot exceed one. Thus 0 < C M < 1 (3.5) Second, again for quasi-incompressibility assumption, the total amount of mass transported into any layer (represented by elements with the same destination index; i.e., rows of the T M ) minus the total amount of mass transported out off that layer (represented by the same source index; i.e., columns of the T M ) has to sum to zero, to conserve the amount of mass in the layer. Since transilient coefficients represent mass fractions, those two facts translate into mass and state conservation conditions (Stull, 1984) requiring that each row and each column of every T M sum up to 1. Hence N ^ C ^ t , At) = 1 for each i (3.6) N J2c{i>j](t,At) = 1 for each j (3.7) for a column of uniform air-mass grid cells. For have the form (Stull, 1993): JV i = i a non-uniform grid, the conservation conditions = 1 for each i (3.8) = 1 for each j (3.9) where [m] is the total mass within any grid cell. In mathematical/statistical terms, mass and state conservation result in T M s that are doubly-stochastic. Finally, dependence of the T M on integration time-step must reflect the real properties of mixing. When At is long compared with the characteristic time-scales of turbulent motions, then the T M approaches a matrix with all elements equal to each other within any turbulent subdomain; namely a well-mixed layer. On the other hand, as A t becomes shorter, all off-diagonal elements approach zero and the diagonal elements approach one; i.e., the T M approaches an identity matrix [/]. Figure 3.4 illustrates those trends using the T M from Figure 3.3 as an example. Assuming that a turbulent sub-domain spans Nsa- layers, the limiting behavior for 22 1 0 i 0 0 | 0 | .5 .5 0 0 0 ! o 1 1 0 0 j 0 | .5 .5 0 0 0 | 0 o i 1 0 ! 0 i "<-*'• [C](At) 0 0 1 0 0 1 o 0 0 1 | 0 i 0 0 0 .5 .5 ! o 0 | 0 0 1 0 0 0 .5 .5 Figure 3.4: Limits of the transilient matrix, [C](At) from Figure 3.3, as A t approaches 0 or oo. corresponding part of the T M , [ C 7 ^ ] , can be written in the following form for a uniform grid: C f c | J ' > « > - { ! IH]} i f A ' - ° (3 ,0 Use of the transilient matrix to classify turbulence The T M can reflect classical idealized properties of turbulence such as homogeneity and isotropy of vertical turbulent motions, as well as stationarity of the turbulence (Stull, 1984, 1993). For these special classes of turbulence, the T M has the following properties: • homogeneity: Cjjj] = C [ i + m J + m ] for all m (3.11) • isotropy: C[v+m] = C[ M -m] for a 1 1 m (3-12) • stationarity: [C](ti, At) = [C](t2, At ) for ti jL ti (3.13) Stationary of turbulence requires more attention since it allows a special composition of T M for multiple or variable A t . For example a matrix for longer time-step can be represented as an inner product of matrices for shorter time-steps: [C](t,At! + A t 2 ) = [C](t,At!) • [C](t + A t i , A i 2 ) (3.14) In general the dependence of the T M on the length of the integration time-step, for the case of stationary turbulence, can be expressed as [C](t,aAt) = { [C] ( t ,A t ) f (3.15) where a is a factor that can be either an integer or a real number. Equations (3.14) and (3.15) should be used with caution if the T M is applied during a period of active external forcing, such as including sources and sinks for passive tracers, or surface forcing for dynamically active quantities like virtual potential temperature or wind speed. Stull (1993) discusses further factors that violate Equation (3.15), including convective structure memory. 23 The early parameterizations of the T M satisfied the so-called exchange hypothesis. Namely, turbulent transport of fluid from layer j to i is to be accompanied by an equivalent transport of fluid in the opposite direction. The original exchange hypothesis is expressed in the following form (Stull, 1984); namely, a symmetric matrix C M = C M (3.16) and reflects the condition of turbulent-transport as imposed by some of the parameterizations of T M . For the case of a modeling framework with non-uniform grid (Stull, 1993), the exchange hypothesis is expressed in the following form: m[J]C[J,i] = m[i]C[i,j] (3-17) where [m] represents the mass of the fluid for each layer. Note, that even though the T M for the non-uniform grid may visually look asymmetric, this asymmetry is superficial and due to the modeling framework, and thus does not represent asymmetry of turbulent transport. 3.1.2 Turbulence statist ics based on the transi l ient mat r i x The T M can not only be used to analyze the classical properties of turbulence, but it can also serve to calculate new types of statistics of turbulence: • Kinematic nonlocal turbulent flux (Stull, 1991b) : Az k N FX[k] = ^ J2 JI (c\iA 3W - CM (3-18) i=lj=k+l is a nonlocal equivalent of classical local turbulent flux at level k for a given quantity S. • Transport spectrum (Stull, 1991b) A k N •= I \ Z ' T >, V \ TS{k,m] = At 1^ 6m,\i-J\(CM E[i] - C[i,j] s b] ) ( 3 - 1 9 ) i=lj=k+l describes the contribution of eddies of sizes mAz to the total nonlocal turbulent flux. • Air transport spectrum (Ebert et al., 1989) Az k N ^ m ] = ^ E ZE 5m,\i-j\(C\Jd ~ C M ) ( 3 - 2 0 ) i=l j=k+l expresses the amount of air movement across level k contributed by eddies of sizes mAz, even though the total air-mass flux must be equal to zero due to mass conservation. • Process spectra (Ebert et al., 1989) k N PS[ktm](At) = J2 5 m , K - j | ( C ^ ] + C [ i , i ] ) ( 3 - 2 1 ) i=lj=k+l provide spectral decomposition, regardless of the associated flux, of total mixing at level k into eddies of sizes mAz. 24 • Mixing intensity (Ebert et al., 1989) k N MI[k](At) = J2 £ ( C y , i ] + C M ) i=ij=k+i describes total mixing intensity at level k for all eddy sizes. Mixing length (Ebert et al., 1989) N 1 ML[k](At) = AzJ2 ^(C[i, fc] + C M ) | i - k\ i=i (3.22) (3.23) describes the average distance that air is relocated due to turbulent motions from level k. Ebert et al. (1989) and Stull (1991a,b, 1993) provide a very detailed and extensive statistical analysis of Transilient Matrices obtained from both observations and numerical experiments. The T M statistics presented above apply to the special case of uniform grid spacing only. Analogous statistics can be derived for non-uniform vertical grid-spacing. In addition to generic transilient statistics presented above, Bagliani (1997a,b), introduced a set of specific statistical metrics, or indices, allowing one to quantify different aspects of asym-metry of Transilient Matrices. First, the author introduces an invariant measure of a state of turbulence represented by T M that is independent of its symmetry. It has the following mass-dependent form for a non-uniform grid: N-l N /([C],[m]) = £ J2 U-i)m[{lCM i=l j=i+l (3.24) or the following dimensionless form obtained by dividing the above equation by the average mass of the grid cells (m) j N-l N J m ([C], [m]) = = Y Y (•?' - »H-]C[U] (3-25) i=l j=i+l where m 1 N -Y (3.26) i=i This invariant, which has the same value for both the upper and the lower triangles of the T M , serves as a basis to calculate a set of four other metrics: • two 'indices of absolute spectral asymmetry' that are always positive, independent of the sign of asymmetry: A 1([C],[m]) = -m/([C],[m]) N-l N YI Y H A J ] - ^fu,^ i=l j=i+l A2([C],[m)) = J([C],[m]) N-l N Y Y im[<] c[i,j]-mb-] cb-.<]iij-' i=l j=i+l (3.27) (3.28) 25 two 'indices of spectral asymmetry' negative or positive depending on the sign of asymmetry: "jV-1 JV J V - l JV A3([C],[m]) = I([C],H) i = l j=i+l j=l i=j+l A4([C},[m}) = I([C],[m]) J V - l JV J V - l JV J2 m [ i i c [ i j ] i j - * i 2 - 2 J2 m w c [ M i b - ^ i 2=1 j = 2 + l j = l i=J + l (3.29) (3.30) A l l of the above metrics have a common feature; i.e., they are equal to zero if T M is symmetric. However, they react differently to the distribution of asymmetry in a T M . Indices A l and A3 are not sensitive to the range of motions represented by transilient coefficients, because they treat all elements equally. Indices A2 and AA give more weight to larger-range transport. 3.2 Parameterizations of the Transilient M a t r i x As introduced earlier, the T M can be parameterized as either an a priori or responsive turbulence closure. Stull (1984) suggested a few different a priori approaches to constructing T M s . Those parameterizations have limited applicability, mostly to the engineering or idealized flows, and as mentioned in Chapter 2, they are not capable of accounting for the wide range of forcing scenarios existing in natural flows. Hence, further discussion here focusses on responsive parameterizations. The application of responsive T3 to atmospheric flow (Stull and Driedonks, 1987; Raymond and Stull, 1990; Stull and Bartnicki, 2000) consists of the following two elements of the turbulence generation-dissipation process: destabilization and stabilization. Destabilization is a buildup of instabilities in the flow by all the various non-turbulence processes; e.g., the potential for turbulence to form due to shear in the mean flow, or due to convective instabilities resulting from external forcing factors such as solar heating of the ground or radiative cooling of cloud top. Stabilization is a response of the flow, mainly by turbulent advection/mixing processes, acting to reduce the instabilities and to lower the turbulent potential accumulated in the flow. The T M represents only the stabilization stage. The destabilization portion of each time-step is crucial to the proper behavior of the T M , because of its nonlocal-in-time character. In discrete/gridded numerical models, variables in the bottom grid are strongly affected by fluxes imposed at the bottom boundary, and are thus sensitive to the thickness of that grid. The same relates to any forced layer such as the cloud top, but it is especially visible in the bottom one, since usually it has the smallest thickness in typical models. (Similar reasoning applies to ocean models, except that the top layer might be considered to be a forced layer.) For example, in a modeled C B L the heat capacity of the forced layer depends not only on the specific heat of air (Cp), but also on its depth. Therefore the relative temperature change of this layer, as compared to any other layer of C B L , depends on its thickness. This, in turn, may lead to unacceptable errors in determining nonlocal gradients. For such reasons, the destabilization process can be split into two parts (Stull and Bartnicki, 2000). The duration of the first part is chosen to be proportional to the height of the bottom forcing layer. The proportionality factor, as determined from numerical simulations, equals 10. The shallower the forcing layer, the shorter is the destabilization stage. The remaining part of the boundary condition is applied after T M is computed, but before T M is applied to do the mixing via Equation 3.1. This split of the time-step mitigates the effect of unrealistically strong gradients on the calculated mixing potential. 26 Two types of responsive parameterizations of T M were previously formulated and applied to atmospheric and oceanic flows. One approach is based on the bulk Richardson number, and the other on a diagnostic form of the T K E equation. 3.2.1 Transilient matrix based on Richardson number Stull (1984) developed a parameterization for the transilient coefficients based on the bulk Richardson number (Equation (1.8)), henceforth called i?#-TM, where the nonlocal equivalent of the bulk Richardson number is g At'iQy A'h r ^ ' i ~ A*jev[(AW)* + ( A W ) * ] T O R * ^ 3 { 6 - 6 L ) The normal bulk Richardson number indicates the dynamic instability of the flow; namely, its potential for becoming turbulent. The nonlocal extension of this Richardson number concept indicates nonlocal dynamic instability and the potential for nonlocal turbulence. This is then used to specify a nonlocal forcing function % i ] = ( ~ ) wM + 1 - w[itj] for i < j (3.32) Larger values of X indicate a greater nonlocal stability, which implies that greater turbulent mixing is needed to undo the instability. This nonlocal forcing factor depends on the turbulence termination value, RT, and on a weighting factor w\hi] = i U° ^ 1 f o r * ^  5 (3-33) { ' \j — i\Az Uo in above equation is a characteristic scaling velocity of turbulence and can be modeled using scaling velocities appropriate for the state of the simulated boundary layer. The values of wui^ approach zero for short time-steps and large separation distances, but are not bounded for large time-steps and short distances. Hence, [w] are truncated to fall below one 0 < w[hj] < 1 (3.34) in order to prevent excessive mixing at large separation distances for longer time-steps. Similarly, the value of [X] is truncated to fall between zero and one 0 < X M < 1 (3.35) in order to prevent mixing at r^j > RT; namely, when the flow ceases to be dynamically unstable. On the other hand, the forcing function is set to one also whenever > RC> where RC is a critical Richardson number, to satisfy turbulence domain limitations. Hence, RC is used to determine the domain of turbulence, and RT is used to terminate the turbulence within the domain. In addition to the above conditions, X[jj] is not allowed to decrease as \j'• — i\ decreases for any i, in order to assure that the strongest mixing will occur between nearest layers. Finally, the differences between transilient coefficients are defined in terms of the forcing functions Cm ~ U ~ = X M for i < j (3.36) k=l which together with mass conservation condition (Equation (3.6)) and the symmetry condition (Equation (3.16)) form a set of A^(A^ — l ) /2 equations that can be solved for a complete set of transilient coefficients necessary to construct a T M . 27 A structurally simpler, although conceptually identical, version of above parameterization was applied by Stull and Kraus (1987) in an ocean model. Apart from modifications necessary to ac-count for specifics of the ocean modeling environment, two changes were made to the formulation of the parameterization. The weighting factor is defined as UpLpAt . . ™M = (j _ ^ 2 A z fort^j (3.37) where LQ is a length scale depending on the type of instability; i.e., it is either the depth of the convective layer for free convection, or the depth of the fluid for wind-forced turbulence. Also, Equation (3.36) is replaced by a simpler one JV Y C M ~ JClhJ] = X[ij] f o r * < 0 (3-38) k = j + l where %;] = ( i - ^ ) (3.3?) Again, one must solve a set of linear equations to obtain the transilient coefficients. A significant change to the way the transilient coefficients are calculated was introduced by Zhang and Stull (1992). They started with the weighting function of Equation (3.33) and the forcing matrix of Equation (3.39). However, instead of being defined by a set of implicit equations, the off-diagonal transilient coefficients are made directly dependent on the forcing function C M = — X M for i ± j (3.40) where is the number of layers in the turbulent domain. The diagonal elements are then calculated from the mass conservation condition (Equation (3.6)); i.e., JV C M = 1 - £ CM (3.41) Thus, the parameterization does not require solving a set of equations for transilient coefficients, and is computationally more efficient than former versions. In summary, all the parameterizations based on the Richardson number are symmetric with respect to indices, when applied to a uniform grid. They depend on a set of three parameters, Uo, R C , and RT, and on a set of turbulence detection conditions depending on those parameters. Zhang and Stull (1992) reported that the best agreement between modeled and observed turbulent mixing in their case study was achieved when Uo = 0.5[m/s], RC = 1.5, and RT = 2.0. However, those particular numbers might not be general enough to apply in the wide range of atmospheric conditions. In particular, the Richardson numbers are significantly different from those used by Stull (1984) and from typical reported values (Stull, 1991b). 3.2.2 Transilient matrix based on diagnostic TKE The other responsive parameterization (Stull and Driedonks, 1987; Raymond and Stull, 1990; Inclan et al., 1996) of mixing potential is based on the turbulence kinetic energy ( T K E ) budget equation, henceforth called the diagnostic T K E - T M . This particular equation is used because it better describes the state of turbulence due to shear instability and convective instability than 28 the Richardson-number method, and it includes viscous damping (dissipation) of turbulence. Computation of T M consists of two major steps: defining a mixing potential (MP) based on the T K E , and normalizing the M P to ensure its conservation properties. In essence, the T M is a normalized M P . Stull and Driedonks (1987) start the derivation of the M P by using the T K E tendency equa-tion, with neglected mean advection and turbulent transport terms, in the form: 9 dtE = -wudzU - wvdzV + -^-w8v - e (3.42) where E is the T K E . Then they define a dimensionless nonlocal analogue of the T K E equation by assuming that instability between two points will generate T K E on the same scale as the distance between those points to remove the instability. B y decomposing the T K E into contributions from eddies of different scales, and by generating a dimensionless form of the equation using the respective T K E , they obtain the following dimensionless equation of T K E in finite-difference form: ^AtEi,j = ( M i j M _ H ^ A W g W ~ _ SiA Eij \ Eij A^z Eij A^z ^ Sv[i] Eu Eitj J ^ 6 ) where AAt denotes a difference over the integration time-step, Eij is the portion of T K E of scale \j — i\Az, and indices i, j represent the properties of the eddies extending from source grid j to destination grid i. The above equation is used to approximate the relative scaling of transilient coefficients, assuming that T K E can be linearly decomposed by spatial scale. The dimensionless form of the T K E equation contains four unknown higher-statistical-order terms: u and v turbulent momentum-fluxes, turbulent heat-flux, and T K E dissipation rate. Three scaling parameters are introduced to parameterize those unknown quantities: a turbulence time-scale To, a critical Richardson number Rc, and a dissipation factor D. Then, the flux terms in Equation (3.43) are approximated using nonlocal gradients of the related variables - ( 3 ,4 , Eij u A'iz (wv)i,• m AW - - r o ^ 7 7 - (3.45) Ei i u A>3z (w8v)hj _ T0A>iev Eitj Rc A>iz and the dissipation term is defined as g j j = D Eij T0 (3.46) (3.47) The right-hand side of Equation (3.43) reflects the contributions of major physical processes that drive turbulent mixing, and Stull and Driedonks (1987) postulate that it can be used directly to define off-diagonal elements of a mixing potential (MP) matrix [Y]. Replacing the unknown terms in Equation (3.43) with their parameterizations yields _ r 0A* (A'V)2 + ( A W ) 2 - - ^ - ^ ( A ^ ) ( A ^ ) To for i ^ j (3.48) with the following suggested values of scaling parameters: T0 = 100[s], Rc = 0.21, and D = 1. 29 The M P matrix defined by Equation (3.48) must be conditioned before it can be used to obtain the T M . First, it is assumed that if the loss terms exceed the production term, then there is no turbulence at all: Y M = m a x ( 0 , Y [ y ] ) (3.49) Second, assuming that large turbulent eddies will stir up a spectrum of smaller scales of turbulence in the regions across which they move (i.e., turbulence cascade), it is required that elements of the matrix are not decreasing when marching from the corners of the matrix toward the diagonal: I m a x ( Y M , Y [ i j + 1 ] , Y [ i + l j ] ) if j < i K 1 Finally, the diagonal elements of M P are defined. Those elements, if within turbulence do-mains, represent sub-grid turbulence. They are defined to be larger than their neighbors in the same matrix row, again due to the turbulence cascade process: Y [ M ] = m a x ( Y [ i ] i _ 1 ] , Y [ M + 1 ] ) + Y r e / (3.51) where the additional parameter, called a reference potential, Y r e y , accounts for internal mixing within layer i. The ratio of the off-diagonal to the diagonal elements has large control on the amount of turbulence and mixing, and a value of Yref = 1000 is suggested, based on sensitivity studies. The M P matrix prepared according to above recipe is ready to be normalized to obtain the T M . A n LQO norm of the M P is defined as I I Y l l o ^ max ( E Y M ) ( 3 - 5 2 ) Then the off-diagonal transilient coefficients are defined from: iii^j (3.53) The diagonal elements of T M are then calculated by a formula identical to Equation (3.41) N — 1 -[i,i = 1 - E C M ( 3- 5 4) and thus the process of constructing T M is complete. In summary, Stull and Driedonks (1987) proposed a responsive and symmetric parameteriza-tion of T M for a uniform grid based on a diagnostic form of the T K E equation, using four scaling parameters. Raymond and Stull (1990) introduced minor changes to the transition from M P to T M in order to apply T3 in a mesoscale weather forecast model. The modifications are done to allow for a non-uniform grid. The LQO norm of of Equation (3.52) is modified to include layer mass: lN \ ||Y||oo = max X l m b ' ] Y M I (3.55) The off-diagonal transilient coefficients based on this new norm have the form ^ " T m u 1 ( 3 - 5 6 ) 30 and the diagonal elements are still defined by Equation (3.54). More significant changes were presented by Inclan et al. (1996) to model turbulent transport across a forest canopy. The modifications improve dependence of M P on integration time-step, and ensure the expected theoretical behavior of T M in the limit of large and small A t . First, Equation (3.48) is replaced by one that better captures T K E associated with nonlocal differences in shear and convective instabilities Y[hj] = (A^U)2 + (AW)2 - ^ ^ - ( A ^ e „ ) ( A ^ ) - D(Ai'h)2 for i ± j (3.57) where the dissipation parameter is redefined to have units of [ s - 2 ] . The suggested values of parameters for above equation are Rc = 0.5 and D = 1 0 ~ 6 s - 2 . Again, conditioning of the M P is defined by Equations (3.49) and (3.50), and is applied first. Next, a renormalization procedure is performed on off-diagonal elements of M P to mimic the behavior of the T M for large and small A t (Equation (3.10)) Y [ « 1 = W T T i f ^ J ( 3 ' 5 8 ) where a parameter (3 = 5 • 1 0 - 4 was determined from L E S simulations. This allows the transilient matrix to approach the identity matrix for small time-steps, and to approach the well-mixed matrix for large time-steps. The diagonal elements of M P are then set to a reference potential Y M = Y r e / (3.59) where Y # e j = 0.5 was determined from numerical experiments. Additionally, for the case of a non-uniform grid, the following weighting procedure, described by Stull (1993), is applied Y M = Y M — 3 - 6 0 The procedure sketched above yields a symmetric mixing potential. Even though the normal-ization can be done according to a simple and fast procedure based on M P norm for symmetric matrices (Equations (3.52), (3.53), and (3.54)), Stull (1994a) finds this approach unsatisfactory because it does not guarantee the proper limiting behavior of T M for large and small time-steps. Instead, Stull (1994a) and Inclan et al. (1996) propose another algorithm, suited for a case of either symmetric or asymmetric matrix. The relaxation of M P towards the T M is done via an iterative procedure suggested by Peterson (Stull, 1995). The scheme redistributes, in a mathematical diffusion-like way, the elements of the matrix, until it reaches the doubly-stochastic form within some error tolerance. It conserves the main features of the M P and, to a certain extent, the relative properties of the matrix elements. The procedure involves a few steps: 1. set the initial guess for the T M to equal the M P : [C] = [Y] 2. calculate vectors consisting of sums of elements in corresponding rows, [R], and sums of elements in corresponding columns, [C): N % = X C M r - mra CIJ] - Z ^ M u i ; i= i 31 3. calculate errors expressing departure of row/column sums from conservation conditions: i=i 3 = 1 4. if both errrow and errco1 are smaller than some prescribed maximum value, e r r m a x , then quit the iteration with [C] representing desired T M , otherwise continue to the next step, 5. calculate a matrix, [ C ] r o w , consisting of [C] with rows independently normalized by R[i] to sum to 1: rnrow fcij] -"-[*] 6. calculate a matrix, [C] c o J , consisting of [C] with columns independently normalized by C[i] to sum to 1: (T'COI 7. replace [C] with an average of [C]row and [ C ] c o L [ c ] _ [ c ] — + [cri 8. go back to step 2. The speed of convergence of the relaxation procedure depends on both the degree of imbalance and number of elements in M P . The value of maximum error depends on desired accuracy of conserva-tion conditions; e.g., Inclan et al. (1996) suggested errmax = 10~ 4. Despite added computations compared to direct calculation of T M , the advantage of this approach is its independence from symmetry of M P . 3.3 S u m m a r y Discrete T M models presented above parameterize nonlocal turbulent advection in physical space as a function of nonlocal mean properties of the flow in physical space, and explicitly simulate nonlocal turbulent advection. Despite the differences in formulations, those models share the following common features: • They are first-order statistical closures, simulating turbulent advection of the mean-state values. • They are responsive parameterizations, where the amount of mixing is related to the amount of instability in the flow; i.e., the models simulate the action of turbulent advection only for conditions where there is a potential for nonlocal/local mixing. • They simulate the combined action of both the smaller chaotic turbulent eddies and the larger coherent structures, resulting in a nonlocal turbulent transport of the properties of the flow between distant locations. 32 • They are not based on the turbulent diffusion concept, but describe turbulent dispersion as an advective-like process. However, if no large eddies exist in the flow, then T M models mimic the small-eddy transport of K-theory. • They are based on limited set of three (for iJg-TM) or four (for diagnostic T K E - T M ) scaling parameters that are determined through numerical sensitivity experiments; namely, ./V2 transilient coefficients are determined from only a small set of pre-defined parameters. • They satisfy the exchange hypothesis (if asymmetry resulting from non-uniform grid is neglected) and represent symmetric transport between model layers; i.e., the amount of mass transported upward between two layers equals the amount of return flow transported downward. • They still assume that eddies have the ability to contribute to nonlocal transport in the entire turbulence domain no matter how short the integration time step of the model is, al-though the measures are taken to limit intensity of long-range transport for short integration time steps. However, the transport range does not depend on the time step. • They are diagnostic, based on the flow conditions at the start of each time step, and do not take into account the history of turbulence; i.e., they discard the information from previous time steps. The research and modeling efforts involving T3 to date, many of them cited at the beginning of this chapter, provided a thorough overview of the behavior of TM-models in different atmospheric conditions. Based on those efforts, some concerns were expressed and suggestions for future development were proposed. 3.3.1 The motivation for this research The exchange hypothesis (i.e., the imposition of matrix symmetry) is considered to be the most questionable constraint of the existing parameterizations. The asymmetry of turbulent trans-port is especially important in the case of convection-driven turbulent domains, as discussed in Chapter 1, which clearly violates the exchange hypothesis. Hence, the need for an asymmetric formulation of the T M has been expressed in the literature. At the time when the first diagnostic T K E - T M model was formulated (Stull and Driedonks, 1987), Stull (1988) presented differences in idealized turbulent dispersion of symmetric and asymmetric T M s , and indicated the need for an asymmetric formulation of T M . Further insight into the asymmetry issue was later revealed by L E S experiments of Ebert et al. (1989), provide estimations of asymmetric T M in convective conditions. Asymmetric T M s were further confirmed in the water-tank experiment of Cenedese and Querzoli (1997). These simulation results can now be used to investigate asymmetry properties of T M s . Stull (1991a) compared the LES-estimated vs. parameterized transilient matrices, and pointed out that the exchange hypothesis was a major deficiency of T M parameterizations. One possible approach to parameterizing asymmetric T M s was proposed by Stull (1995), who suggested using convective available potential energy ( C A P E , Equation (4.27)) in place of the convective-forcing term in Equation (3.57). Bagliani (1997a) followed this direction and encoun-tered a problem with convergence of the Peterson relaxation scheme, which required artificial compensation of elements in the M P matrix in order to ensure convergence, and which did not admit much further development of the model. The goal of developing a stable and reliable asymmetric T M parameterization is a challenge that stimulated this dissertation research. 33 Other minor deficiencies of the existing parameterizations are reported in different studies. The most comprehensive analysis, done by Stull (1991a), indicates an unrealistic vertical distri-bution of mixing intensity in the C B L ; namely, an exponential decay with height and very weak mixing in the entrainment zone near the top of the C B L . Raymond and Stull (1990) also find the intensity of mixing unrealistically high near the bottom of turbulent domain. This characteristic of T3 leads to super-adiabatic profiles extending far into the upper part of C B L , as reported by Zhang and Stull (1992); Cuxart et al. (1994). Interestingly, Zhang and Stull (1992) report that i ? B - T M gives more realistic profiles in the surface and mixed layers of the C B L than the diagnostic T K E - T M . Another deficiency indicated by Stull (1991a) is related to the balance be-tween off-diagonal and diagonal elements of the T M , depending on renormalization of off-diagonal elements and the pre-defined (but somewhat arbitrary) value of the reference potential Yref. In summary, the weaknesses in the current T3 parameterization method are: • The imposition of symmetry in transilient matrix. • Unrealistic vertical distribution of mixing intensity in the convective boundary layer. • A d hoc parameterization of diagonal elements in the transilient matrix. • Infinite transport-range (a common feature of nonlocal closures, as reviewed in Chapter 2). This is the motivation for the current research. The introduction of higher-statistical-order methods into T3 nonlocal closure is explored in this dissertation, in order to investigate whether asymmetric T M parameterizations are possible, and whether the other deficiencies of former T3 parameterizations can be mitigated. In the next chapter, this will be implemented via a prognostic parameterization for the T M . To support the new parameterization, the following concepts are introduced: • Turbulent-transport mean-flow state operators: a transilient tendency-matrix as a tendency operator, and the transilient matrix as an integral operator. • A turbulent-advection model based on a virtual turbulent-transport eddy concept. • Integral measures for mean-flow instabilities: convective potential-energy and shear potential-energy. • Symmetric and asymmetric components of turbulence-forcing by mean-flow instabilities. • Two alternative models for turbulence forcing by mean-flow instabilities. • A prognostic model for the nonlocal turbulence-kinetic-energy budget. • A prognostic model for the 3rd statistical moment of the vertical velocity. In subsequent chapters, the parameter space of this candidate closure method is explored, and resulting turbulence and mean-state simulations are compared against results from other models. 34 Chapter 4 A prognostic model for the transilient matrix For any conserved tracer E, if we let [£](£) be the vector of concentrations at discrete grid points in a column, then T3 describes the evolution of tracer with time by [S](t + At) = [C](t,At)[2](t) (4.1) where [C] is the transilient matrix (TM). However, this can be equivalently written in more generalized form as a tendency equation dt[E](t) = [C](t)[E](t) (4.2) where [C] is an operator that will be called here the transilient tendency-matrix ( T T M ) . The T T M is not dimensionless like the T M , but has units of [ s - 1 ] . Although Equation (4.1) de-scribes turbulent transport that is nonlocal in time (according to the terminology described in Appendix A ) , Equation (4.2) describes turbulent transport as local in time. While previous transilient parameterizations focused on describing [C], the approach taken here is to parameterize [C] as a function of a prognostic form of the bulk T K E equation. This is analogous to higher (one-and-a-half) order local closure, which carries a prognostic equation for local T K E . Thus, this chapter will present a higher-statistical-order nonlocal closure. First, the transilient-turbulence-theory generalization is explained using a system of turbulent-transport equations for any mean-flow quantity, which yields the corresponding turbulent-trans-port operators. The introduction of a conceptual virtual turbulent-transport eddy follows. Pre-sented next is an extension to nonlocal static-stability criteria. Then a set of prognostic equations is presented that allows the calculation of eddy-wise virtual turbulent-transport velocities. Finally, a new definition of mixing potential (MP) is proposed together with a new way of transforming the M P into a T M . Throughout this chapter, as in Chapter 3 (see page 18), the finite-difference view of a fluid column represented by N layers is assumed, with a discretized system of conservation equations in both space and time. As before, it is assumed that nonlocal turbulent motions between two layers of the atmosphere (i and j) result from nonlocal mean-flow instabilities between the same layers. As in traditional Transilient Turbulence Theory (T3), closure of this one-dimensional discrete column model is based on only the information available about the flow; namely, the time-dependent mean properties of the flow field, as well as the independent variables of time (£) and space (z). 35 4 . 1 Prognos t ic transil ient model ( P T M ) concepts 4.1.1 Turbulent- transport mean-flow state operators To demonstrate how Equations (4.2) and (4.1) are related, start with the budget equation for S, and focus on the contribution from the vertical-turbulent-transport term: dtZ(z,t) = -dM(z,t) (4.3) where E is the mean tracer concentration, £ is the turbulent fluctuation from the mean, and w£ is turbulent kinematic flux. Integrate it over the depth Az^ of any grid cell: Ot^[i](t) = ^ — (4.4) Now, decompose the turbulent flux at each interface between two grid cells into source-oriented contributions from turbulent transport between all layers in the fluid: k N N k = Z Z * M s [ 7 l ( * ) ~ Z Z ^ M S b ) W (45) j = l i = k + l j = k + l i = l where W[ij] is some conceptual turbulent-transport velocity from layer j to i. The first term on right-hand-side corresponds to upward turbulent flux, the second to the downward. There is no turbulent flux at the boundaries, hence the following conditions complete the definition in Equation (4.5): ^(*)[0]=^(<V]=0 (4.6) Next, replace the turbulent fluxes in Equation (4.4) with Equation (4.5), by defining w£(t)t0p = w£(t)[i] and w^(t)i,ottojn = (£)[,_!], to obtain after simplifications the following budget equation: j ^ i k+i (4.7) where first term in the brackets on the right-hand-side represents the incoming turbulent-flux, and the second term represents the outgoing. Finally, consider the system of all Equations (4.7) for all grid cells in the column. The result, in abstract matrix form, is dt[E](t) = [C}(t,.. .)[E](t) (4.8) where [C] is a generic mean-flow state-tendency operator, defined as follows: Cp „(*,. .) = ( (4.9) It is a 2nd order tensor (matrix) explicitly acting upon vector [E]. In general [C] depends on other variables and parameters in addition to t, so the ellipsis symbol is used here for brevity. If [C] is quasi-stationary during time interval A t , than Equation (4.8) can be explicitly solved: [E}(t + At) = e® At[E}(t) (4.10) 36 But, the exponential of the operator in Equation (4.10) acts like an operator itself. Define this new operator as [C](t ,At , . . . ) = e[C](t--)A< Thus, Equation (4.10) can be abbreviated as [~](t + A i ) = [C](t ,Ai, . . .)[~](i) (4.12) However, one can recognize Equation (4.12) as identical to the basic T3 framework Equa-tion (4.1), and Equation (4.8) as being identical to the new P T M framework Equation (4.2). The key is Equation (4.11) which relates the transilient tendencjf-matrix (TTM) to the T M . The above results mark a distinction between former version of T3, and the new parameteri-zation of the T M proposed here. As presented in Section 3.2.2, both the M P and then the T M were directly parameterized as a function of an integration interval At. Below, a parameterization for a new M P (which only roughly corresponds to the former M P , as it describes a potential for turbulent transport) is proposed that is semi-independent of A t (except for the limited turbulent-transport range, as will be explained later in this chapter). That new M P is then used to define the T T M ([C](t)). Equation (4.11) is used next to obtain a At-dependent T M ([C](t, At)) . There are two important consequences of this new approach: 1) implicit state conservation, which is built into Equation (4.2), eliminates need for a separate state-conservation Equation (3.9), and 2) it will be shown later in this Chapter that there is no need to parameterize the diagonal elements of the M P , as was done in former T3 parameterizations. 4.1.2 Virtual turbulent-transport eddy concept As described in Chapter 1, the turbulent transport of any quantity in the atmosphere is realized by cumulative action of eddies of different scales, and with different levels of coherency. Such transport is a complicated function of u(x,y, z,t), v(x,y, z,t), and w(x,y, z,t). Since T3 de-scribes net nonlocal turbulent-transport for only the vertical component of a turbulent field, only w(x,y, z,t) is considered. The concept of a virtual turbulent-transport eddy ( V T T E ) is intro-duced here. It allows one to approximate the turbulent-transport of a quantity using a simple bimodal model of vertical turbulent-transport velocity. Consider the vertical turbulent-transport of a single conceptual eddy that is exchanging air-mass between two horizontal layers of the atmosphere. Those layers are separated by vertical distance h, and spread over a fixed horizontal area A, such as the area of a single grid-cell in a N W P model. Thus, the volume spanned is B = hA. Volume averaged 2nd, W2, and 3rd, W3, statistical moments can be defined as follows: W2(B,t) = ^ JJJw2(x,y,z,t)dB (4.13) B W3(S, t ) = ^JJJw3(x,y,z,t)dB (4.14) B where dB is a dummy variable of integration. One may separate the volume integrals into two parts; namely, those with upward motions in sub-volume Bup(t) and those with downward motions in sub-volume BdOWn{t): W2(B,t) = i JJJw2(x,y,z,t) dB + ^ JJJ w2(x,y,z,t) dB (4.15) W3(B,t) = ^JJJw3(x,y,z,t)dB + ^JJJw3(x,y,z)t)dB (4.16) BuP(t) Bd0wn(t) 37 Now, assume that one can approximate the turbulent-transport velocity in sub-volume Bup(t) by some characteristic or average velocity Wf(B,t), and in sub-volume Bd0Wn(t) by another char-acteristic velocity u>i(B,t). Replacement of w(x,y,z,t) in the above integrals with the corre-sponding characteristic velocities yields: W2(B,t) = wi{B,t) + ^ ^ wt2(B,t) (4.17) W3(B,t) = ^4^w,3(B,t) - B d 0 " n ( * W ( / J , t ) (4.18) JD JD Note the sign convention introduced above; namely, W[(B,t) represents magnitude of upward velocity, w±(B,t) represents magnitude of downward velocity, and both are always positive. B it*) Now define a partial volume of upward motions by P|-(.B,£) = u £ ' and a partial volume of downward motions by Pi(B,t) = Bd°™n^, Then assume that the air mass inside of the box is conserved; i.e., mass-flux out balances mass-flux into the box. This yields a set of four equations: 1 = P^B^ + P ^ t ) (4.19) 0 = Pi(B,t)wi(B,t)-Pi(B,t)wl(B,t) (4.20) W2(B,t) = Pi(B,t)wi2(B,t) + Pl(B,t)wi2(B,t) (4.21) W3(B,t) = Pi(B,t)w]3(B,t)-Pl(B,t)wi3(B,t) (4.22) The above set of equations closely resembles the simple bimodal top-down/bottom-up turbulent advection model proposed first by Wyngaard (1987) and then applied to nonlocal closure by Zilitinkevich et al. (1999). The main difference is that, as presented here, the averaging is done over a volume, not over a horizontal plane. Note that if one knows the value of W2(B,t) and W3(B,t), then the Equations (4.19) through (4.22) form a closed set and can be solved for P T ( B , i ) , Pi(B,t), wi(B,t), and Wi(B,t): W,t) = \ - , a W 3 < B ' * > (4.23) ' 2 2v/W32(B,t) + 4 W 2 3 ( £ , t ) V 1 Pl{B,t) - U « E L ^ (4.24) 1 2 2y/W32{B,t)+4W2*(B,t) MB,t) = W3(B-t)+c^ ;r4W— (4-25) w i ( B i i ) _ - m ^ ^ ^ t y . g M ( 4 2 6 ) The above V T T E equations define an approximation of the turbulent-transport field based on four quantities: characteristic virtual turbulent-transport partial-volume ( V T T P ) of updrafts and downdrafts, Pj(B,t) and P | ( 5 , t ) ; characteristic virtual turbulent-transport velocities ( V T T V ) of updrafts and downdrafts, IOJ(£?,£) and wi(B,t); and their relation to basic properties of the vertical turbulent-velocity field (2nd and 3rd statistical moments as approximated by W 2 ( 5 , t) and W 3 ( £ , i ) ) . Equations (4.23) through (4.26) can be analyzed to show that the mathematical properties represented in those equations correspond to the anticipated physical properties of the real tur-bulent vertical-velocity field; namely, (1) a finite value of W2 implies finite T K E associated with vertical turbulent transport, and (2) a finite value of W3 introduces a measure of asymmetry into the upward and downward motions. Specifically, positive W3 results in narrower but more rapid updrafts (as observed in convective boundary layers), while negative W3 indicates the opposite 38 effect. W3 = 0 results in symmetric motions as expected for neutral and stable layers, where the asymmetry is either negligible or not observed. A real drawback of these equations is their dependence on the 3rd statistical moment of vertical turbulent-velocity. Usually an estimate of that quantity requires at least a 2nd-order closure. Nevertheless an attempt is presented in Section 4.2 to introduce a simple model for W2 and W3 that is sound enough to simulate desired properties of turbulent-transport field, but without complexity of full 2nd (or higher) order turbulence-closure. 4.1.3 Extension of nonlocal static stability criteria Stull (1991c) remarked on a potential inaccuracy of his nonlocal static-stability scheme, presented on page 6, due to neglecting the overshoot effect. By neglecting overshoot one ignores the fact that an air-parcel's inertia could allow for nonbuoyant penetration into stable regions. That may lead to underestimating the depth of turbulent layers. He suggested that a more sophisticated nonlocal definition of static stability could be utilized to correct this effect. One can extend Stuff's nonlocal static stability criteria by borrowing the concept used rou-tinely in detecting convective instability for thunderstorm forecasting. The concept is based on the convective available potential energy ( C A P E ) . The form of C A P E suitable for turbulence prediction is the one using virtual potential temperature (e.g. see Stull (1991b)) fZLOC a CAPE(zLFC, zLOC) = / 7Jnl®v(z) ~ e*(z)] dz (4.27) JZLFC V \Z) where ©^ is virtual potential temperature in the cloud, ©^ is virtual potential temperature in the surrounding environment, Z^FC is the level of free convection at which cloud air starts to have positive buoyancy with respect to its environment, and z^oc is the limit of convection (called also a level of neutral buoyancy (LNB) (Bohren and Albrecht, 1998)) where cloud air stops having positive buoyancy with respect to its environment. According to observations of convective clouds, at least part of the cloudy rising air, the updraft core, is relatively undiluted, with little mixing with its environment. If one then uses the environmental virtual potential temperature at the L F C as an estimate of ©^ in the updraft core, then the definition of C A P E (Equation (4.27)) can be rewritten as fZLOC a CAPE(zLFC, zLOC) = / T ^ T T [ Q V ( Z L F C ) - ®v(z)] dz (4.28) In reality the physics of convective clouds is much more complex, due to the mixing processes at the top and the sides of forming clouds, water phase changes inside of the cloud, and wind shear in the cloud environment. Nevertheless, the C A P E concept still serves as an indicator for the vertical range of convective motions in convective clouds and thunderstorms. Even though C A P E is most commonly used to determine potential for convective updrafts-, the properties of the C A P E integral allow it to be used equally well for the opposite case; i.e., for cold air sinking. A n example is when a relatively colder layer aloft is created in the atmosphere, e.g. due to radiative cloud-top cooling. In such a case, the C A P E integral is applied in reverse direction for cold air sinking. The change in sign of integration dummy variable dz together with the negative [@V(ZLFC) ~ ®v(z)]> due ^° * n e negative buoyancy force, results in a positive sign of C A P E for convective downdrafts as well. Hence, C A P E positiveness is characteristic for any convection-driven motions, independent of whether the convection is created due to warm air rising or cold air sinking. Under closer inspection one can notice a correspondence between C A P E and the nonlocal static stability definition of Stull. First, the C A P E integral starts at the L F C , corresponding to 39 the local maximum in virtual potential temperature profile, from where air parcels are lifted to detect nonlocally unstable layers. Second, the C A P E integral stops at the L O C , corresponding to the level of neutral buoyancy in detection of unstable layers. Finally, the C A P E integral correlates the sign of buoyancy force with the direction of considered motion, and the nonlocal static stability detection uses the sign of the buoyancy force in qualitatively the same way. In fact, Stull's definition of nonlocally unstable layers is equivalent to finding the layers between the points having positive C A P E with respect to each other. The C A P E stability concept can be further expanded by taking into account the inertial overshoot effect. The parcel moving under the force of buoyancy wil l convert its potential energy into kinetic energy, and in the extreme case when there are no energy losses due to mixing or air drag, its kinetic energy at the L O C will be equal to the C A P E . When the parcel overshoots, the reverse energy conversion takes place, namely its kinetic energy is converted into potential energy until the overshooting air parcel looses its inertia. In the extreme case, when no other energy losses take place, the motion will continue to the level at which the negative counterpart of the C A P E integral in the deceleration zone balances the positive C A P E . In order to quantify this effect one can generalize the C A P E definition (Equation (4.28)) in the following way CAPE*(d, s) = max ( I 7£-r[Qv(s) - Sv(z)] dz , o] (4.29) \Js ®v(z) J where s corresponds to initial source level of integration, and d corresponds to final destination level of integration. In Equation (4.29), d and s can be any two levels within the grid column. Using the generalized C A P E definition, one can identify nonlocally unstable layers, including the overshoot effect, as those where for each point in the layer there exists a pair of points, bottom (zb) and top (zt), such that either CAPE*(zb, zt) > 0 or CAPE*(zt,zb) > 0. Figure 4.1 shows the results of applying both Stuff's and the generalized C A P E algorithm to the detection of convectively unstable layers in the atmosphere. The algorithm based on generalized C A P E detects a larger extent for the unstable layer due to the included overshoot effect, while the unstable layer according to Stuff's method is entirely included (as a sub-domain within) in the generalized-CAPE detected unstable layer. While Stuff's method under-estimates the depth of convectively unstable layer, the generalized C A P E method over-estimates its depth, because it does not account for energy losses (e.g., the energy cascade to smaller eddies that is associated with turbulent drag) in the conversion between kinetic and potential energy. So far in this section, only the buoyancy force was considered as a source of vertical motions. However, the kinetic energy of vertical turbulent motions can be provided also by shear instability. If the source of kinetic energy is other than the buoyancy force, and if the buoyancy force is solely working as an inhibitor of vertical motions, one could supplement the generalized C A P E with an approach similar to convective inhibition (Bohren and Albrecht, 1998), used for convective systems, and define simply the convective potential energy (CPE) CPE(d,s)= fd -J—[ev(s)-ev(z)]dz (4.30) Js Vv[Z) The C P E can be interpreted as follows: buoyancy will enhance parcel's motion in the direction that decreases its C P E ; will oppose the motions that increase parcel's C P E ; and will not affect the motions that do not change parcel's C P E . In other words, C A P E indicates how much energy is available before the air parcel rises, while C P E indicates how much potential energy was consumed after it rose. Using the above definition, one can provide a quantitative nonlocal static stability condition equivalent to Stuff's method (page 6), but including the overshoot effect. Using the same require-40 z 120 100 80 60 Stull 40 C A P E 20 273.5 274 274.5 275 275.5 Figure 4.1: A n illustration of nonlocal stability detection using Stull's vs. the generalized-CAPE algorithm. A warmer layer is introduced into sub-adiabatic virtual potential temperature profile. Vertical lines represent the depth of the convection-driven turbulent layer as detected using: Stull's nonlocal static stability method (limited by stars); and the generalized-CAPE method (limited by boxes). These are the statically unstable regions. ments as in Stull's method, the C P E nonlocal static stability can be determined in the following sequence of steps: 1. unstable: if for each point in the layer there exists a pair of points, bottom (zb) and top (zt), such that either CPE(zb, zt) > 0 or CPE(zt, zb) > 0, 2. stable: layers that are not unstable, and such that for each pair of points (t, 6) in the layer both CPE(t, b) < 0 and CPE(b, t) < 0, 3. neutral: layers that are not unstable, and such that for each pair of points (t, b) in the layer CPE(t,b) = 0, 4. unknown: any stable or neutral layers if parts of the profile are missing and therefore C P E is undetermined. The C P E nonlocal-static-stability criteria not only provides the means to determine the static stability of atmospheric layers, but it also can be incorporated into a new definition of the Richardson number to determine the dynamic stability when shear production is also present. To do that, a concept can again be borrowed from convective-systems forecasting, where a specialized form of the bulk Richardson number (BRN) (Moncrieff and Green, 1972; Weisman and Klemp, where U is defined as the difference between density-weighted mean wind speed taken over the lowest 6km of the atmosphere and an average wind speed taken over the lowest 500m. The B R N is used to anticipate the type of thunderstorm (e.g., a multicell storm is likely for BRN > 30, and supercell storms can grow for 10 < BRN < 40), but its original definition and criteria cannot be readily applied to turbulence. 1982) BRN = CAPE 0.5t72 (4.31) 41 The B R N concept can be adapted to turbulence. The numerator of the bulk Richardson number (Equation (1.8)) is really a crude approximation of C P E , and it is approximately equal to —2 CPE if a linear virtual potential temperature profile is assumed. The Richardson number reflects the ratio of buoyancy to shear forcing terms in the T K E Equation (3.42), and the specifi-cation of those forcing terms defines the particular form of the Richardson number. Both forcing terms originate from the tendency equation for the velocity perturbation, and the term under the C P E integral is equivalent to a vertically-integrated buoyancy force. One could apply the same approach and specify an equivalent vertical integral of shear force - call it the shear potential energy (SPE) - which would have a form f d SPE(d,s) = {-udzU - vdzV} dz (4.32) J s where u and v are fluctuations of the horizontal velocity components. Approximating the per-turbations of velocity by the difference between mean velocity at the starting point and mean velocity at the environment through which the parcel is moving, the above definition can be expressed as SPE(d, s)= f {-[U(s) - U(z)]dzU - [V(s) - V(z)]dzV} dz (4.33) Ja It can be readily integrated to yield SPE(d, s) = \{[U(d) - U(s)f + [V(d) - V(s)}2} (4.34) This S P E is always positive, as expected for shear instability. Also, it is equal to 1/2 of the denominator of bulk Richardson number (Equation (1.8)). Therefore, for the case of an idealized near-linear virtual temperature profile for stable layers, one can approximate the bulk Richardson number using C P E and S P E as: _ -CPE(d,s) =. R b " SPE(d,s) ( 4 3 5 ) The close similarity between expressing the bulk Richardson number in terms of a finite difference approximation of the gradient Richardson number, and expressing it for stable layers using C P E and SPE, suggests that one can define yet another form of Richardson number, call it a P E Richardson number, based entirely on C P E and S P E -CPE(d,s) R p E = SPE(d,s) ( 4 3 6 ) Further, the similarities indicate that one can expect RPE to have properties comparable to the bulk Richardson number: allowing determination of dynamic stability, and possessing limits of critical and terminal values. Finally, one can expect that C P E and S P E can be used to define buoyancy and shear forcing terms in the kinetic turbulence Equation (3.42), as proposed in the next section. 4.1.4 Convers ion between potential and kinetic energy A simple way to model the evolution of turbulence, and the properties of V T T E , is by means of the Turbulent Kinetic Energy ( T K E ) . Before proceeding any further, recall simple facts from elementary mechanics. Consider a body, either solid or not, with constant mass. The force working on the body causes it to accelerate and gain kinetic energy. The acceleration of the parcel is equal to the force divided by the mass of the body; i.e., specific force. The gain of 42 specific kinetic energy will be proportional to the distance the body moved and the acceleration of the body due to the acting force: (specific kinetic-energy gain) = (specific force) * (distance) The rate at which the specific kinetic energy changes is proportional to the gain divided by the time over which the gain of kinetic energy was achieved: . „ , . . . (specific kinetic-energy gain) (specific kinetic-energy gam-rate) = — -(specific force) * (distance) (time) Distance divided by time is equal to velocity, hence the rate at which the body gains kinetic energy, under the influence of constant specific-force, is: (specific kinetic-energy gain-rate) = (specific force) * (velocity) (4.37) Now look at the acceleration process from slightly different point of view. If the body is moved vertically by a force, one can assign to the body a property called potential energy; i.e., the amount of energy that can be transformed into kinetic energy. Once the body starts to move it looses its potential energy. The acceleration, i.e specific force, of the body is proportional to the loss of potential energy divided by the distance over which the loss happened: / -c r ^ (specific potential-energy loss) (specific force) = (4.38) The last two equations can be combined to yield a simple, yet useful, formula describing trans-formation rate of potential-energy into kinetic energy: / - f l , • ,. • . x (specific potential-energy loss) (specific kmetic-energy gain-rate) = — r * (velocity) (4.39) (distance J Equations (4.37) through (4.39) may serve as templates that describe conversion rate by which some form of 'turbulence' potential energy characteristic for some forcing mechanism in the flow field is transformed into turbulent kinetic energy. This concept will be used in the next subsections to help to parameterize the TKE effects on VTTE. However, in addition to TKE, the definition of the VTTE also requires a model for the 3rd statistical moment of turbulent velocity, w3(x, y, z, t). The simple physics considerations above are not valid for this term. However, one can utilize an analogy between TKE and 3rd statistical moment of vertical velocity to complete the VTTE definition, as described next. 4.2 Prognos t ic equations for W2 and W3 In this section, simple parameterizations for W2(B,t) and W3(B,t) are presented that allow calculation of VTTV's and VTTP's from equations (4.23) through (4.26). Also, the nonlocal turbulent transport of any quantity in the volume B (between fluid layers separated by distance h) will be approximated. Simulating the evolution of turbulence, and especially the asymmetry of the vertical turbulent-transport field of the VTTE, is the main goal of this new approach. In order to achieve that goal, the equation for the 3rd statistical moment of vertical turbulent velocity is introduced, and 43 the equation for T K E is utilized in its prognostic form. As was done in the former T M model, the equations are based on nonlocal properties of the mean flow, and result in nonlocal turbulent transport of the mean-flow properties. To develop the parameterization, one starts by formulating the physical properties of the proposed model. As presented in Chapter 1, generation of turbulent motions depends on the existence of turbulence-production mechanisms, which in turn depend on the mean-flow state of the fluid. The turbulence itself can also support generation of turbulent motions through the 'cascade' mechanism of transfer of energy through the spectrum of turbulent motions; i.e., in gross simplification, the turbulent motions at larger scales contribute to the generation of turbulence at smaller scales. On the other hand, turbulence is also self-limiting and dissipative. Namely, once it exists, it acts to remove the source of instability that originally generated the turbulence. It also acts to eventually eliminate turbulence altogether, again mostly through the turbulence energy cascade to the smallest (mm size) eddies, where it is ultimately dissipated at the molecular level. The state of turbulence depends on the balance between production and dissipation; dom-ination of one of them will lead to either intensification or reduction of turbulence intensity. Therefore any model of turbulence should consist of at least: turbulence production by mean flow instability, turbulence energy-gain due to energy transferred from larger scales (energy gain from the cascade process), and turbulence energy-loss due to energy transferred to smaller scales (energy loss to the cascade process). This can be expressed in generic form: dtT = (T-forcing) + (T-cascadeGain) + (T-cascadeLoss) (4.40) where T is T K E , (T-forcing) represents production by mean-flow instabilities, (T-cascadeGain) represents T K E gain from the cascade process, and (T-cascadeLoss) represents T K E loss to the smaller eddies. The full T K E tendency equation (see Equation (2.8)) has additional terms expressing turbu-lent transport and pressure redistribution of T K E , which are not considered here. A n attempt was made in this work to parameterize those effects as a common term (in a way equivalent to local closures). However, the parameterization of that term proved to only marginally improve the prediction in statically neutral and stable conditions, but was a source of considerable nu-merical instability in convective conditions. Therefore it is omitted in this dissertation, and left for future work. Earlier T M parameterizations also neglected these terms (see Chapter 3). The asymmetry of vertical turbulent motions, present in convective turbulent flows, is usually considered only in higher-order closures. For example, in 2nd-order closures it is usually ap-proximated via constant skewness, which admits parameterization directly as a function of T K E . Here a prognostic model based on buoyancy forcing and decay of W3 is proposed instead. This approach, in generic form, is dtS = (S-forcing) + (S-decay) (4.41) where S is the 3rd statistical moment of the vertical-velocity field, (S-forcing) represents buoyancy forcing, and (S-decay) represents decay. In the remainder of this chapter, the focus will be on the asymmetry of vertical motions. Two velocity scales are used in this work. The first one is a T K E velocity scale ([715]), which is independent of any anisotropy in the flow because the T K E is based on a sum over all three Cartesian directions: (4.42) The factor 2 comes from the fact that T K E in any one direction is equal half of the variance in that direction, and the factor 3 comes from partitioning T K E equally between the 3 Cartesian 44 directions. Although larger [US] indicates larger T K E , the [715] contains no information about asymmetry of the turbulent motions in any one direction. The second scale, defined here as an asymmetry velocity scale ( [ .45] ) , expresses the magnitude of asymmetry of the turbulent motions in the vertical direction: •AS[i,3]W = V S M ( T ) ( 4 - 4 3 ) Namely, this scale is a measure of skewness of the turbulent vertical-velocity distribution. The transport properties of turbulent motions can be related to the turbulence. As illustrated in Chapter 1, turbulent transport is advective in nature and turbulent motions carry fluid (and all its properties) from one location to another. Therefore turbulent-transport scales and properties are closely related to the scales and properties of turbulent motions. Under that assumption, one can relate the turbulent-transport scales to properties of V T T E , as above. Furthermore, it is assumed that turbulent motion at a given scale results in turbulent transport at the same scale. Hence the following relations hold W2(B,t) = j J T M ( i ) (4.44) W3(B,t) = SM(t) (4.45) where the factor 2 appears as before. Substituting these into equations (4.23) and (4.26) yields: ^[.jlW - 2 2Jsfij]W + i T 3 j ] W (4.46) 2Jsy*) + § T 3 j ] W (4.47) fT[<j](*) These P and w values can now be used to build the prognostic T3 parameterization. 4.2.1 Turbulence generation by mean-flow instabilities Shear and convective instabilities in the mean flow generate turbulence for as long as turbulence itself does not remove the cause of the mean-flow instability. Hence, one can define a nonlocal specific-force, Fjjj] (£), acting at scale — on the air due to instabilities between layers j and i. It is a sum of shear F^^(t) and buoyancy FB^(t) specific-forces (expressed in units of [m/s2]): ^ M W ^ M W + ^ I W ( 4 - 4 8 ) The shear specific-force represents turbulence generation when the mean-flow velocity varies with height. A n air parcel moving from level j to i that maintains its horizontal momentum will 45 have a velocity different from its new surroundings, which results in a turbulence kinetic energy. As proposed in the previous section, one may say that the shear-related potential energy of a parcel at level j with respect to level i is equal to S P E defined by Equation (4.34), which in the discrete modeling framework has the form s p E U t j _ ( A ^ I O W A W f f ( 4 4 9 ) Since the change takes place over a distance Al<h, applying the template Equation (4.38) to S P E yields: Q SPEHi](t) F M w = n^ Si (4-50) The shear specific-force is symmetric with respect to indices i and j, and it is always positive. Namely, it results in production of T K E . The buoyancy contribution to the specific-force term represents the process of turbulence production/loss whenever the virtual potential temperature varies with height. Even though initially Stull (1995) proposed C A P E to estimate asymmetric mixing potential, the formula for C P E (Equation (4.30)) is used here as an expression of the potential energy associated with motions from layer j to i. The discrete version of the formula is: CPE[hj](t) = ZT^TTT {©«[>](*) - ®v[k](t)} Az[k] sgn(i - j) (4.51) k = j K y v [ k ] \ z ) where additional function, sgn, is defined as follows - {r1 Ixtl <4-52> It is used here to preserve the sign of the original integral. Again the change takes place over a distance A%,Jz when parcel is moving at vertical velocity TVS. Applying template Equation (4.38) yields: where a form of the critical Richardson number RC>T is introduced here as a parameter that allows one to change the balance between the shear and buoyancy forcings. i?C,T includes the hysteresis effect (Equation (1.9)) and is equal to: 1) the critical Richardson number, RC, if eddy is at least partly outside of the turbulent domain, or 2) the terminal Richardson number, RT, if eddy is entirely inside the turbulent domain. Buoyancy specific-force is source-oriented, and therefore asymmetric with respect to indices i and j. It is solely responsible for introducing asymmetry into vertical turbulent-transport. It is either a production or loss term in the T K E equation depending on the type of nonlocal stability of the atmosphere between layers j and i, as discussed in section 4.1.3. So far, the specific forces acting on a parcel moving unidirectionally from level j to i were considered. However one may recall from Chapter 1 that turbulence can have a vortical character, and parcels moving from one level to another will be accompanied by a return flow in the opposite direction, in order to conserve mass. Therefore one might consider the forces acting in the eddy circulation together, how the balance of those forces contributes to the production/loss of T K E , and how it affects asymmetry of the motions. Two alternative models are considered next: a VTTE-wise forcing model and a parcel-wise forcing model. 46 V T T E - w i s e model One may apply the V T T E concept to the forces acting on the circulation between layers j and i, henceforth called VTTE-wise model. First, define average specific forces acting in the sub-volumes of the upward motions Ff ^ ^ (£) and downward motions F^ ^ ^ (t). Then recall that the equations for T and S are originally obtained by multiplying the perturbation equations by the appropriate powers of perturbation velocities. Using w-\ and W[ as surrogates of perturbation velocities for the upward and downward sub-columns, then integrating the forcing terms over a V T T E yields the following total forcing terms (the source/destination indices and time-dependence are temporarily dropped for clarity): (T-forcing) ~ P T w T F T + P ^ J F I . . (S-forcing) ~ 3(PTV^T " A«>l2*i) where factor 3 is retained from the original equations for S. Now, define the following 'new' forces: 1) a 'symmetric' force acting on the circulation as Fs = F t + F j (4.55) and 2) an 'asymmetric' force acting on a circulation as FA = F T - F x (4.56) where in both cases half of each matrix is redundant. The last two equations can be inverted to obtain F T = ( F 5 + F^ ) /2 p l = ( F 5 _ F ^ ) / 2 C4-57) Introducing the above forces might seem artificial at first. In order to show the physical properties of those forces, define F j and F j in terms of mean-flow specific-forces expressed by Equation (4.48): v - tit] <«•> One can show, by replacing the forces on the right-hand side with Equations (4.50) and (4.53), that the complete form of those forces for mean-flow instabilities is F 5 r ,+s _ (A^[/( t ) ) 2 + ( A W ( t ) ) 2 gcAUeu(t) V z A a w Qualitatively, the symmetric force is equivalent to former T K E - T M model for the M P (Equa-tion (3.48) or (3.57)) if the summation term divided by the nonlocal height difference is taken as a measure of inverse M:^QV. The asymmetric force is a function of curvature of the Qv profile, and is positive in typical convective P B L and negative in typical top-cooled cloud layers. Both forces are symmetric with respect to indices, and although half of each matrix is redundant, it simplifies notation if these forms are preserved in further considerations. Equations (4.57) together with V T T E Equations (4.46) and (4.47) can be inserted into Equa-tions (4.54) to obtain the following tendency equations: (T-forcing) ~ TVS (S-forcing) ~ T 2 T 3 /2 S 2 + 16T3 F 5 ^ 3S p S (4.61) 47 The above equations proved to be computationally very sensitive to small values of T and S. Hence, the next simplification is done by defining a skewness parameter as _ sgn(§)S W = ( |T)3 /2 (4.62) and using it to eliminate T and S from inside the brackets in Equations (4.61): (T-forcing) ~ 1\S (S-forcing) ~ T YA _|_ sgn(S) Sw (4.63) where Sw is a parameter that must be calibrated against canonical case studies. Parcel-wise model According to Equations (4.63), T K E is not generated in free-convective layers for those eddies equal and larger than the depth of any statically unstable layer, defined in terms of Stull's nonlocal stability definition, because F 5 = 0, and hence it does not allow for convective penetration. This deficiency can be remedied if one considers a parcel-wise model. That is equivalent to either dropping partial volumes in Equations (4.54) or using the template Equation (4.37) separately for the upward and downward sub-volumes of V T T E : (T-forcing) (S-forcing) (4.64) Repeating the former derivation procedure for this new forcing leads to the following tendency equations (T-forcing) (S-forcing) TiS V 8 1 S 2 + 9 6 T 3 F 5 , 3 V 2 s ' V A 4 x / 6 T 3 / 2 T ^ 4 T 3 / 2 E T ( 2 + | g ) F ^ +3 V 8 1 S ; + 9 6 T 3 F 5 (4.65) 8 T 3 -that after substitution using Sw definition have a form (T-forcing) (S-forcing) TVS V / 4 + 5 l u 2 j p 5 sgn(S) Sw ^A T (2 + SW2)¥A + sgn(S)Swy/^ + Sw2Fs (4.66) Those equations have distinctly different form from Equations (4.63). The T tendency depends now on F" 4 , and this model can potentially allow for the inertial overshoot effect in penetrative convection, in the sense of the generalized C A P E definition of static stability introduced in Section 4.1.3. Summary The VTTE-wise and parcel-wise models contain the following features: • S generation is possible only if T K E T > 0. • In both models the T forcing can be expressed in the form of template Equation (4.37) with TVS as the velocity scale. 48 In both models the S forcing can be expressed in the form analogous to template Equa-tion (4.37) with T as the velocity-variance scale. The VTTE-wise model of T forcing does not depend on asymmetric force. In consequence, as already mentioned above, the eddies deeper than the depth of statically unstable lay-ers (defined in terms of Stull's nonlocal static-stability) are not generated, and convective penetration through a capping inversion is not simulated. The parcel-wise model allows for convective penetration through asymmetric forcing in the T tendency-equation, and increas-ing the value of Sw parameter will lead to increased penetration depth. However, the depth of penetration is limited by the condition that (T-forcing) > 0; i.e., Fs > —FAsgn(s)|S™ ; where the fraction is bound in an interval (—1,1). Increasing Sw has opposite effects on T K E forcing between the two models. In the V T T E -wise model, the forcing monotonically decreases down to a zero limit, as Sw increases. In the parcel-wise model the forcing behavior depends on the combination of signs of Fs, FA, and S, and it need not be bound as Sw increases. However, it can be shown that for an eddy-size corresponding to the depth of the free-convective layer (in terms of Stull's nonlocal definition, where CPE^j] = —CPE^^), T forcing is equal to SWF^^ = SvRC^Jif • Hence, it is reasonable to assume that Rc should define an upper limit for Sw, if physical values of CPE are to be considered. In both models, in the limit of Sw —> 0, the T forcing is exclusively dependent on Fs, and the S forcing is exclusively dependent on FA. VTTE-wise model and par eel-wise model have different limits as Sw —* 0. The magnitude of the VTTE-wise forcing is half that of the parcel-wise forcing. In order to match those two models in limit of zero Sw, the forcing in VTTE-wise model will be multiplied by 2. (This choice is arbitrary, but it also matches T forcing models to the forcing of the former T K E - T M model, as mentioned in the discussion after Equation (4.60)). S tendency may be non-zero in any turbulent layer. However, it is commonly accepted that non-zero skewness is an important property of convective layers only, and it is usually not even considered in stable and neutral conditions in typical closures up to 2.5-order. Therefore § forcing and Sw will be set to zero outside of convective layers; i.e., when buoyancy forcing has a damping effect on turbulent motions. There is no compelling evidence of significant local skewness reversal inside of convective layers, and their skewness depends on the character of the dominant forcing. However, both models of S forcing allow for significant local skewness reversal, specially at the edges of the convective layers when Fs is negative and strong compared to FA. On the other hand, negligible skewness in stable and neutral layers suggests that it is not important when buoyancy has a damping effect on turbulence. Therefore, it will be further assumed here that the S forcing is considered only if the T forcing is positive, with S forcing set to zero otherwise. Both forcing models depend on sgn(S), which in general does not have to be the same as the sign of FA forcing. In practice, it might lead to undesired numerical effects. It is safe to assume that the sign of S is the same as the sign of FA, which is a reasonable assumption except for very rapid transitions between stable and convective regimes that rarely appear in the real atmosphere. 49 Taking into account the above considerations, the T K E specific-force, Fj^j^t), becomes: • VTTE-wise model: * T M ( ' ) = * T (4.67) • Parcel-wise: Fni,j](l) = KT V / 4 + £ (4.68) where Kj is the T K E specific-force parameter. The asymmetry specific-force, Fg[jj](i), becomes: • VTTE-wise model: FS[i,j](t) = { 0 if FJ[id] > 0 i f F T [ l j ] < 0 (4.69) • Parcel-wise: ^S[ij](*) = (2 + 5,1(j2)F-4[jJ-](t) + Sw^/l~Vs^sgn(¥\j](t))Fs[i^](t)\ if F T [ y ] > 0 if FT[iJ] < 0 (4.70) where the Ki is the asymmetry specific-force parameter. The forcing terms in the tendency equations have the form: (T-forcing) = TVSM FT[i>j] (t) (S-forcing) { i J ] = T[jj] F s [ i ) i ] (t) This completes the derivation of mean-flow instability forcing terms in the tendency equations for T and S. These two models, due to their different behavior under increasing 5^ and the common limit for Sw —> 0, can be used to more thoroughly investigate the influence of skewness on the evolution of turbulence and properties of turbulent transport. (4.71) 4.2.2 Turbulence energy-cascade So far, forces acting on the circulation between levels j and i due to properties of the mean-flow field at (or between those levels) were formulated. Next, consider how large turbulent eddies generate smaller scales of turbulence. This generation mechanism corresponds to the turbulence energy-cascade, as explained in Chapter 1. T K E gain from the cascade process Shear production by the larger eddies is usually the major mechanism responsible for smaller-scale turbulence generation, and the only information that is available so far in this model is in T[jj](t) and S[ij]. Hence, those two quantities are utilized by analogy to shear production by the mean flow. To simulate the down-scale turbulence-energy-cascade, a simple mother-daughter model is considered, where mother is a larger eddy feeding energy to smaller daughter eddies. 50 Consider a mother V T T E extending between layers I and J . B y analogy to Equation (4.49), sub-scale shear potential-energy (SPE) for the mother V T T E can be estimated using its updraft and downdraft velocities as K [ M + ^ [ M ) 2 ( 4 7 2 ) which after applying Equations (4.47) yields 9 s[i,J\ (4.73) 8 T M ] Next, consider a daughter V T T E extending between layers j and i, entirely enclosed between layers I and J , and gaining energy from its mother due to the turbulence energy-cascade. First, one needs to estimate the sub-scale S P E associated with the daughter. The numerator of Equa-tion (4.72) has a form that may by considered as a velocity structure-function (mentioned in Section 1.1.2 of Chapter 1) associated with the mother V T T E . Therefore, one can assume that the sub-scale S P E associated with daughter should: 1) be a fraction of that associated with the mother, 2) decrease to 0 as the daughter size decreases to 0, and 3) approach the mother's S P E as the daughter's size increases. To approximate that behavior, a simple normalized linear function is utilized: A^z A^h (4.74) which is bound to the interval (0,1). Then, a daughter's specific-force due to the sub-scale shear can be written, using template Equation (4.38), as 9 § [ M , 4 + 3 T M 8T, VA A'h A^z \A<JZ\ 9 § [ / , J ] 8 X AJ'Jz (4.75) Next, the conservation of T K E in the energy-cascade must be considered: 1. Let A K / ' J ) = \I — J\ be the number of resolved daughter scales that are smaller than the scale of mother [I,J]. The linear function in Equation (4.74) is multiplied by 2 / iV( / , J ) to sum up to 1 over all scales. 2. Let N^.J}ij be the number of daughters in the vertical direction with a scale having an index difference of \i — i.e., having the same size if resolution nonuniformity is neglected. The normalization factor is assuming that each daughter at specified scale draws the same share of energy from the mother. 3. Next, consider the number of daughters in the horizontal direction. The area covered by a daughter is a (A'iz/A1'^)2 fraction of the area covered by the mother; hence, the daughter's integrated production rate should be weighted accordingly. Then, the complete normalization factor to Equation (4.75) (which mimics turbulence-energy conservation during cascade) has the form: Ai'h N ( i , j ) N ( i , j ) \ £ i , j z i(j-i) (4.76) and the complete specific-force contribution to daughter [i, j] from mother [I, J] due to the sub-scale shear has the final form {A<iz) 2 r ' 9S [ f | J ] ( t ) | 8 (4.77) 51 where KjG is a T K E cascade-gain parameter. In order to complete the procedure, and fully use the template Equation (4.39), Fji^t) must be multiplied by the appropriate velocity. When the T K E velocity scale is used with the term describing conversion of potential energy into kinetic energy, the result is: (T-cascadeGain)^ j } = TVS[itj](t)FG{iJ](t) (4.78) The above equation approximates the turbulence energy-cascade in terms of the T K E generation due to shear production by larger eddies. T K E loss to the cascade process A commonly used approximation of turbulence dissipation based on a single dissipation scale like the Taylor microscale (e.g., Tennekes and Lumley, 1972) or as used in local closures (e.g., Equa-tion (2.12)) is not applicable for a multiple-turbulence-scale problem considered here. Instead, a spectral-like approach is implemented (Tennekes and Lumley, 1972); namely, it is assumed that the T K E loss to the cascade process is proportional to some power of T K E and inversely proportional to the length scale at which the loss takes place. In order to preserve the proper dimensions, the term can be defined as (T-cascadeLoss)^-] = - K ° L ^ 0 ^ (4.79) where KjL is a T K E cascade-loss parameter. 4.2.3 Asymmetry decay Asymmetry decay can be defined as ( s . d e c a y ) „ , = _ y f ^ M ( y ) r * ( 4 8 0 ) where K? is an asymmetry-decay parameter. The sgn function (Equation (4.52)) is used to preserve the sign of asymmetry, and make sure that the decay is always acting against asymmetry. 4.2.4 Summary After combining all terms together, the resulting equations for T^^t) and §[jj](£) tendency (with US replaced using Equation (4.42)) are: = Tf ( T M ( * ) ) 1 / 2 + *JJ5](*)} - K$L ( T l A ^ r (4'81) % ] W - % j ] W *S [ i j ] ( * )--K s | (4.82) These are the desired equations, that can be used to predict evolution of T K E and asymmetry of the circulation eddies extending between layers j and i. In summary, two prognostic models for the T K E , [T](i), and the 3rd statistical moment, [§](*), of turbulent circulation were derived above, and the virtual turbulent-transport partial-volumes and virtual turbulent-transport velocities were approximated. 52 4.3 M i x i n g potent ia l Following the approach used in the former T M parameterization, the concept of mixing potential (MP) is used as an intermediate step towards obtaining the T M . The intention here is to define the M P matrix as a prototype of the transilient tendency-matrix, [C], from Equation (4.2). Therefore, the new M P matrix has a different interpretation than in the former T3 parameterizations. Here, it contains information about the potential for turbulent transport, resulting from the turbulence response to flow instabilities. 4.3.1 An MP parameterization The multi-step procedure to define M P starts by utilizing V T T V ' s (defined by equations (4.47)) in its initial form [Yo](t): Y ..(*) = / [ i i > j ( u p w a r d m o t i o n s ) ( 4 g 3 ) °[«J] | w^ jf i < j (downward motions) 1. Each element of the T M matrix describes the fraction of air arriving at a destination grid-cell from a source, due to the cumulative action of all eddies at the corresponding scale, while Yo[ij](t) describes only single V T T E . Hence, the M P must include this information. If the horizontal area covered by the V T T E extending between layers j and i is proportional to (Al'Jz)2, then the factor representing the number of those eddies over a fixed horizontal area should be proportional to l/(A%i3z)2. Two functions can then be built for each layer k that counts only those fractions with nonzero elements of Yo[k,i](z): Fr^[k}{1) f ° r updrafts (sources below and destinations above), and Fr^(t) for downdrafts (sources above and destinations below): Z7V U \ yfc-1 s9n(Y0[k,i](t)) , y^jV s g n ( Y 0 [ i fc)(t)) / A O A \ *rT[fe]W = 2-i=l ( A * ^ — + l^i=k+l—(S^p— ( 4 - 8 4 ) EV f+\ sgn(Y 0 | M 1(i)) s g n ( Y 0 [ i k ] ( t ) ) , . *rl[k]{t) = L,i=k+i ( A ^ ' -Z)2— + 2-,i=i —r^f2— (4.85) where N is the number of vertical layers in the model, and the sgn function is defined by Equation (4.52). Now, the renormalized elements of YFr[k^(t) can be written in a source-oriented form: Y ••(*) = / (A4P Frl(t)Yo[i,j}(t) if * > 3 (upward motions) F r [ h j ] \ Frm(t)Yoii,j}(1) i f * < J (downward motions) It can be shown that the factors in front of Y 0[ij](t) sum to 1 over all elements for the selected layer; hence, they form a normalized distribution-density function. 2. The matrix of normalized V T T V ' s is transformed into a mass-flux matrix by multiplying each element by air density at the source YMF[i,j](t) = YFr[i:j](t)pb] (4.87) 3. The range of transport is limited by assuming that any source elements (from height-index j) with V T T V ' s too weak to reach their destination (height-index i) should be discarded: YLR, ,(t) = ! YMF[i,dt) i f Yo[ij](*)At > (A^z - z{i]/2 - zm/2) lJJJ^ |^  o otherwise ^ ' 53 However, the discarded elements should not be lost entirely. The concept, equivalent to the mother-daughter turbulence energy-cascade in Section 4.2.2, is applied here. Equa-tion (4.88) is applied to all [YMJ?] elements marching from the matrix corners towards the main diagonal. Then, the discarded mother element is redistributed between all daughters (except those on the main diagonal, representing unresolved daughter eddies) using the formula: ^ N M - I 1 J P T Y M F M ( 4 8 9 ) Namely, the discarded elements contribute equally to each daughter scale, and at each daughter scale they contribute equally to each daughter. Thus, a new matrix, [ Y L H * ] ( £ ) , is created, with both limited transport-range and compensation of discarded elements in preserved transport-scales. Note, that limiting transport range changes the geometry of [YLR*](t) compared to the [YMF]^), and that the geometry of [ Y L , R * ] ( £ ) depends on A i . That almost completes the procedure of defining a prototype M P matrix that mimics turbulent transport of air mass. However, before [ Y £ , H * ] ( £ ) can be used further, one last step remains; namely, one must ensure mass conservation of the M P . This is one of the major changes between the former T3 and the new parameterization proposed here. The mass conservation Equation (3.8) and state conservation Equation (3.9) (already implicitly satisfied by Equation (4.1)) are replaced by a mass-flux conservation equation: Ym (*) = E Y [ M (*) f o r e a c h k ( 4- 9°) » 3 This equation says that all the mass flowing out of layer k to all destinations i must be equal to the mass arriving at layer k from all sources j. The numerical procedure used to ensure mass-conservation is presented in the next section. The M P , [Y][jj](i), is defined in such a way as to represent the partial flux of mass transported by turbulence from layer j to layer i, and can be used to calculate the mean-state tendency operator to find the desired T M . 4.3.2 Mass- f lux conservation a lgor i thm The following iterative algorithm (a method similar to the Peterson algorithm described on page 31), referenced later as [Y] balancing, is used to relax'the non-conservative M P into one that is mass-flux conservative. Define the sum of column elements, COL^, (left side of Equa-tion (4.90)), the sum of row elements, ROW[k], (right side of the Equation (4.90)), and the average of those two sums, AVG^, all corresponding to a layer k: COL[k] = J2Y[ltk](t) i i ? 0 ^ [ f c ] ^ E Y [ M ( i ) 3 [k] + ROW[k] 2 Next, recalculate column and row elements according to the following formula AVGlk] for each j Ytk ,-i = Ynu , - i — - — -J 1* Jj l*JJ R 0 W COL{k] + ROWl AVG[k] = -for each i Y [ i ) f c ] = Y [ i i f c ] [*] AVG[k] COL[k] 54 One can verify, that after adjustment, the sum of row and column elements for layer k satisfies Equation (4.90). However, the adjustment does not affect the total mass-flux attributed to that layer; i.e., AVGuk] does not change during the adjustment. The adjustment modifies the elements of the matrix, and has a side effect of changing the sums of the other rows and columns associated with the layer k. Therefore, it must be repeated to march through all layer indices, to apply the adjustment iteratively. Two types of iterative procedures may be considered: (1) a simultaneous procedure where all sums are first calculated for all layers and then applied, and (2) a sequential procedure where the adjustment is done immediately after the sums are calculated and therefore adjustment for next layer uses the already recalculated elements. Numerical tests, performed as part of this research, have shown that the sequential algorithm is numerically stable, and computationally more effective and accurate, and therefore was chosen in this work. The iterative algorithm requires a convergence check. Therefore the error, Err, is calculated after each iteration Err = \ROW[k] - COL[k}\ k and the procedure terminates when Err < Err MAX, where Err MAX is a predefined desired computational accuracy. Err MAX influences computational efficiency of the algorithm and must be chosen carefully as a compromise between accuracy (i.e., mass conservation) and computational efficiency. Both considered procedures have a negative side-effect of attenuating the asymmetry of the final matrix compared to the initial one. Different approaches were also investigated, such as a conservation-condition-constrained minimization, in order to preserve as much of the initial asymmetry as possible. However, they were either computationally less efficient (by at least two orders of magnitude of computation time), or did not result in conservative mixing potential. 4.4 Transi l ient ma t r ix The transilient tendency-matrix (TTM) is derived here from the M P . Then, the T T M is trans-formed into the transilient matrix (TM). 4.4.1 Transil ient tendency-matr ix In order to derive the T T M one must first define its diagonal elements. Those elements correspond to outgoing mass-flux so one can define a new intermediate matrix [C^] as follows: Recall that [Y] defines the partial turbulent mass-flux from source to destination layer. Since the state Equation (4.4) is defined in terms of kinematic fluxes, subsequent steps involve transforming [Cd] back into kinematic form. That is achieved by: 1. dividing each element of the matrix by the air-density at the destination C f c [i,j](t) = — Cd[i,j](t) (4-92) 2. dividing each element of [<C&] by the thickness of destination layer c[i,ii(*) = A I - C * M ( ' ) ( 4 - 9 3 ) 55 Since Azjj] • within a grid column of unit cross-sectional area is equal to the mass of the layer, m^, the entire procedure of transforming [Y] into [C] can be summarized by a single equation: f — £ * Y f J L ,-!(*) if i = j C'«M^S> ^ (4-94) This T T M can be used in Equation (4.11) to calculated the T M . The resulting T M satisfies both the state and mass conservation constraints of T3, provided that the M P conserves mass flux. 4.4.2 Transforming [C](t) into [C](t,At) The transilient matrix [C](t, A i ) is calculated using Equation (4.11), in which the transilient tendency-matrix [C](t) is the argument of an exponential function. The simplest way to calculate [C] is to use a series expansion of the exponential function 2 0 0 n ex = l + x + — + ...= ) — 2! ^ n\ n= 0 and compute the further elements of the series until the series converges to the result at a desired accuracy. However, this approach is limited by slow convergence, and should be avoided. It is mentioned here only to show that the solution yields a matrix composed of real numbers, whereas some numerical schemes will have tendency to generate complex numbers. Those, if small, may be discarded as numerical errors, or a better accuracy scheme can be sought. A n eigensystem method is used in this work, that is more computationally efficient than exponential expansion, and is sketched in Appendix D. 4.5 Notes on T 3 statistics A number of different turbulent statistics based on the T M were presented in Section 3.1.2 of Chapter 3. As a consequence of the limited turbulent-transport-range assumption, the statistics that offer eddy-wise analysis (transport spectrum, air-transport spectrum, process spectrum, or mixing length) are much more sensitive to the integration time-step because larger-eddies might be excluded from the T M . Hence, they do not offer significant insight into general properties of the T M and are not used in this work. More meaningful values could be obtained if one uses a T M based on some characteristic time-scales of turbulent motions. For example, one could use the time scale based on the size and TVS of the most energetic eddies, which would offer more valuable information than that based on A i . Some modification to the statistics are required to fully account for nonuniform vertical grid-spacing and density variation. Hence, a new form of Equation (3.18) for the kinematic nonlocal turbulent-flux was derived, and is used further in this work: A N k N LFCJ i= i j=k+\ The analogous correction could be applied to the other T3 statistics if desired. 4.6 S u m m a r y A new approach, the prognostic transilient model (PTM) model, to parameterize the tran-silient matrix was introduced in this chapter, based on higher-than-first-order statistical nonlocal 56 turbulence-closure. The approach employs mean-state equations and implements an asymmetric model for vertical turbulent transport. The asymmetric model parameterizes vertical-turbulent-transport properties of the turbulence field, based on prediction equations for T K E and the 3rd statistical moment of vertical-velocity perturbations. The virtual turbulent-transport eddy ( V T T E ) concept is utilized to transform predicted properties of the turbulence field into char-acteristic vertical-turbulent-transport velocities. Those, in turn, are used to define a mixing potential that serves as a prototype for the transilient tendency-matrix. Finally, the transilient tendency-matrix is converted into the transilient matrix. In summary, the following assumptions were made, in addition to the initial discretized view of the fluid, in order to build the different elements of the T M model: 1. Assumptions related to mean-state operators: • the transilient tendency-matrix is assumed to be sufficiently stationary to calculate the transilient matrix using Equation (4.11). 2. Assumptions related to the concept of a V T T E : • turbulent-transport velocities in volume B can be approximated by characteristic ve-locity scales Wf(B,t) and wi(B,t), in Equations (4.17) and (4.18), • motions in volume B conserve mass, in Equation (4.20). 3. Assumptions related to the prognostic equations for [T](£) and [§](£): • turbulence is isotropic with respect to spatial scales, in Equation (4.42), • horizontal scales of turbulent motions correspond to vertical scales of turbulence forc-ing, in Equation (4.42), • T K E is equally partitioned between its directional components, in Equation (4.42), • the 2nd and 3rd statistical moments of V T T E correspond to the predicted T K E and the 3rd statistical moments, respectively, of the turbulence field, in Equations (4.44) and (4.45), • the skewness parameter has the sign of the asymmetric specific-force in VTTE-wise and parcel-wise forcing models, in Equations (4.67) through (4.70), • asymmetry is considered only for convective layers and only for positive buoyancy forcing, in Equations (4.69) and (4.70), • the daughter shear potential-energy (SPE) is a linear fraction of mother S P E , in Equa-tion (4.74), • all daughters at any one specified scale draw the same shares of energy from the mother, in Equation (4.77), • the horizontal area covered by the V T T E extending between layers j and i is propor-tional to (N'iz)2, in Equation (4.76), • T K E loss to energy-cascade process is proportional to the 3/2 power of T K E , and inversely proportional to the scale of the eddy, in Equation (4.79), • decay of asymmetry is proportional to the 4/3 power of the 3rd statistical moment of vertical turbulent-velocity, and inversely proportional to the scale of the eddy, in Equation (4.80). 4. Assumptions related to the M P : 57 • each V T T E contribution to turbulent transport from any source-height is inversely proportional to the V T T E area, in Equations (4.84), (4.85), and (4.86), • V T T V ' s can be used to determine the range of turbulent transport, in Equation (4.88). P T M parameterization of the T M utilizes the following eight dimensionless parameters, whose values can be determined through studies of canonical cases, sensitivity studies, and calibration procedures that are presented in the next chapter: RC - critical Richardson number, RT - terminal Richardson number, Sw - skewness parameter, K- - T K E forcing parameter, K~][ - asymmetry forcing parameter, K%G - T K E cascade- gain parameter, KjL - T K E cascade-loss parameter, K§ - asymmetry decay parameter. For comparison (e.g. see reviews by Garratt (1992); Sorbjan (1989); Stull (1991b)), one-and-a-half order local statistical closures, which have been well established in the scientific literature for decades, utilize typically at least three parameters mostly related to length scales. Second-order local statistical closures utilize typically from four to eight parameters, depending on complexity. It is not easy to classify P T M model in terms of closure-order. From the point of view of T K E prediction, it falls into the one-and-a-half order category. But asymmetry prediction is typically reserved for at least second-and-a-half order category. Nonetheless, the P T M model considers a spectrum of nonlocal effects that are completely neglected by most local closures. A l l things considered, it seems more justified to treat P T M as one-and-a-half nonlocal statistical closure. The final set of equations used in the P T M must be solved in the order listed in Appendix D. This is the prototype model that will be utilized in the next chapters. The calculation of the T M , as presented in this chapter, is a computationally expensive task. This increased computational expense should be offset by increased accuracy of the nonlocal turbulence parameterization, as was indicated in Table 2.1. However, some computational savings are possible by including additional non-turbulence terms of the mean-variable tendency equations into the transilient tendency-matrix. This is possible as long as those terms are general for all mean-flow quantities and linear with respect to those quantities. Vertical advection or vertical filtering are good candidates, and combining them into the mean-state tendency operator may save computations. 58 Chapter 5 Estimation of P T M parameters The goal of this chapter is to investigate the parameter space for the predictive transilient theory from Chapter 4. Recall that eight dimensionless parameters were used to formulate the prognostic transilient model (PTM) . Wi th the marriage of such radically different closure concepts (statistical and nonlocal) into one theory, it is anticipated that there will be conflicts between the effects of different parameters. For some atmospheric cases, the higher statistical-order parameters might dominate the calibration, at the expense of those nonlocal-closure parameters that might be pushed away from their optimal values. For other atmospheric cases, the nonlocal effects might be better captured in a different part of the parameter space. By mapping this parameter space, one can better understand potential for future efforts in nonlocal higher-order statistical closure theories. First, the rationale and advantages of canonical case-studies are presented. Then, the available data sets are identified. Next, the experimental design is laid out. Finally, sensitivity studies are used to map the parameter space, and to identify the optimum calibrated values. The evaluation of any turbulence closure is challenging. Ideally, one would like to implement the closure into a three-dimensional numerical weather-prediction (NWP) model, and assess its influence by verifying the forecasts against atmospheric measurements. However, such an ap-proach is difficult due to complex, often nonlinear, interactions between different components of the N W P model. Another, more common approach, involves testing turbulence closures in a one-dimensional (column) version of the N W P model with significantly simplified components. That has its limitations as well, since a number of assumptions has to be made about external forcings, and extensive interpolation/smoothing of observational data is required both to initialize the model and to compare against the N W P results. Measuring turbulence in the real atmosphere is also difficult due to the difficulty of measuring three-dimensional properties, the fact that weather is uncontrollable and unrepeatable, and the need for long time-series of measurements in order to create reliable statistics (Lenschow and Stankov, 1986). The advances in numerical modeling of turbulent flows created another source of evaluation 'data'; i.e., large-eddy simulation (LES). LES allows repeatable, controlled experiments of the flow properties and turbulence-generating instabilities through simple external forcing parameters such as surface heating, surface roughness, and geostrophic forcing. Hence, it can limit the number of parameters influencing the comparison with turbulence-closure models to a small set. It also allows one to extract turbulence statistics at any point of the simulation, thereby providing both initial conditions and verification data for turbulence-closure models. L E S has long been utilized to investigate properties of turbulence closures (e.g., Deardorff, 1972; Moeng and Wyngaard, 1989; Andren and Moeng, 1993). Although L E S cannot entirely replace verification with real atmospheric data, it provides a valuable idealized input for turbulence-closure developers. L E S models are not capable of exactly reproducing all of the features of atmospheric flows. 59 Table 5.1: Archetypical P B L cases of LES runs, the output from which are used to compare with P T M and with other P B L models. (The barotropic cases have imposed constant geostrophic wind with height, while the baroclinic cases impose shear in the geostrophic wind. A l l cases are barotropic, unless noted otherwise.) Case Description E K M Idealized Ekman layer. Shear driven. Neutral static stability. No inversion. OOWC Realistic neutral B L . Forced convection. Capped by weak inversion. OOSC Realistic neutral B L . Forced convection. Capped by strong inversion. 24F Realistic unstable B L . Pure free convection. Capped by strong inversion. 05WC Mixed convection (moderate free + forced). Capped by weak inversion. 03SC Mixed convection (weak free + forced). Capped by strong inversion. 05SC Mixed convection (moderate free + forced). Capped by strong inversion. 24SC Mixed convection (strong free + forced). Capped by strong inversion. 15B Baroclinic atmosphere. Moderate free + forced. Strong capping inversion. 24B Baroclinic atmosphere. Strong free + forced. Strong capping inversion. Ebert Realistic unstable B L . Pure (weak) free convection. Capped by weak inversion. The accuracy of those models depends on factors such as resolution, numerical approximations, or sub-grid-scale parameterizations (Pope, 2000). Some of the reported deficiencies of L E S models are: departure from logarithmic wind profile in the neutral surface layer (Andren et al., 1994), or problems in modeling of evolving stable layers. Nevertheless, both L E S models used in this study were verified against observations (e.g., Moeng, 1984; Moeng and Wyngaard, 1988; Schmidt and Schumann, 1989), and showed satisfactory performance. Hence, the results of the L E S simulations are treated in this work as a sufficient approximation of the true atmosphere, and accurate enough to perform numerical experiments with the P T M model. 5.1 Ava i l ab le L E S cases Two sources of L E S data are used: 1) the LES database ( M M M - N C A R , 2003), created at the U S A National Center for Atmospheric Research (NCAR) for planetary boundary-layer (PBL) model evaluation, and 2) the Ebert case (Ebert et al., 1989) providing estimates of the Transilient Matrix (TM) from a different LES model. Details of the published LES database for P B L model evaluation (Ayotte et al., 1996) are given in Section B . l of Appendix B , including the forcing conditions and procedures used by those investigators to generate initial conditions and evaluation data. The Ebert case (Ebert et al., 1989) is described in Section B.2 of Appendix B. Those archetypical cases are briefly summarized in Table 5.1. In brief, the L E S database for P B L model evaluation provides two sets of data for each L E S case (cases E K M through 24B in Table 5.1): 1) a set of initial conditions consisting of initial mean profiles of horizontal wind components U and V, potential temperature O, and passive tracers B and C; and 2) a set of evaluation profiles including mean profiles of U, V, 0 , B, and C , as well as mean profiles of turbulent fluxes for momentum wu and wv, heat w6, and passive tracers wb and wc. The Ebert case (last in Table 5.1) provides initial conditions for potential temperature 0 , and a set of estimates of the transilient matrix. These data sets allow one to run a single-column P B L model and compare the model results against those from the L E S simulations. 60 5.2 P T M experiment design During parameterization development, the prognostic transilient model was implemented into a simple single-column model (called also P T M ) with 7Y layers of variable depth, [Az], and a constant integration time-step, At. The following additional elements were introduced into the P T M to match the L E S forcing conditions, and to provide the extra elements necessary for P T M parameterization: 1. As with the L E S cases, the mean state of the atmosphere is assumed to be horizontally homogeneous, allowing mean horizontal advection to be neglected. 2. As for the L E S cases, the influence of the horizontal pressure gradient on the horizontal mean-wind field is simulated in terms of geostrophic-wind forcing. 3. The surface-layer parameterization of Louis (1979), based on Monin-Obokhov similarity, is used to describe surface fluxes and their influence on the bottom model grid. The scheme was modified to be driven by prescribed surface heat-fluxes. 4. The P T M requires knowledge of the vertical density profile in the atmosphere. The model is assumed to be in hydrostatic equilibrium; i.e., no mean vertical motions are simulated, vertical accelerations within turbulent updrafts and downdrafts are assumed to be negligibly small, and the vertical air-density distribution is set to that of the standard atmosphere (Holton, 1992). The above elements of the P T M provide a minimal set of forcing conditions that approximate conditions in a horizontally statistically-homogeneous atmosphere, and allow for straightforward implementation of external forcings, initial, and boundary conditions from the L E S cases. Addi-tional details of numerical approximations used to solve P T M equations are given in Appendix D. The small integration time-step and fine vertical resolution of the L E S runs are much finer than the values typical for P B L models. On the other hand, one would like to ensure a fair representation of the evolution and resolution of turbulent flow, as well as to provide a common representation base among the different L E S cases. Therefore, the integration time-step, At, and vertical resolution, Az, for the P T M model are defined based on the large-eddy turnover-time, T, and the P B L depth, Zj , given by: 1. At= I T , 2. Az = ±zi. L E S case-specific r and z, are listed in Table B.2 in Appendix B . The resulting P B L discretization falls into the range of values typical for mesoscale N W P models, and allows for fair analysis of P T M behavior. ; In order to match resolutions of the L E S and P B L models, the L E S initial and final profiles were interpolated to the case-specific PBL-model grids, and used for initialization and evaluation. This completes the procedure of transforming results from the L E S database into a form suitable for P B L modeling. The Ebert case requires a different setup, and will be discussed in Section 6.2. Having a limited number of cases to both calibrate and evaluate the P T M model, one must decide how to use them effectively. In the approach taken here, the L E S database is divided into two groups. The first one consists of three cases: E K M , O O W C , and 24F. They share the common property of simulating the limiting types of turbulence generation listed in Table 5.1. Hence, they allow one to isolate and investigate the influence of different parameters of the P T M on its behavior. Those canonical cases are used to calibrate the P T M . The second group, 61 comprising the remaining cases, is used to evaluate the P T M and to compare it with other P B L models. The Ebert case is treated separately from LES database, and it is used specifically to discuss estimated properties of transilient matrices. 5.3 Parameter-space explorat ion procedure The P T M uses a relatively large set of parameters, and it is not feasible to determine all of them at the same time. The approach taken here is to start the calibration from those L E S cases that affect the smallest number of significant P T M parameters, and then use those already-determined parameters as fixed values in more complex cases. Hence, the first case to be used is E K M , involving only shear-generation, dissipation, and the energy cascade of turbulence. Then, the influence of Richardson numbers on turbulence damping in statically stable capping inversion is determined using case OOWC. Finally, case 24F is used to find the parameters influencing the asymmetry properties of strong convective turbulent-transport. The P T M exhibits transient oscillations in mean and turbulent-flux profiles, due to oscillations in the transilient matrix. Hence, the P T M outputs are temporally averaged (similar to the averaging applied to L E S results as described in Section B . l of Appendix B) . The simulation time of the P T M model is extended for another r /2 , and the final profiles are averaged over T. Up to 10 different variables (winds, temperature, fluxes, etc.), depending on the L E S case, are predicted in each simulation. As in many models, the optimal parameters are found as a compro-mise between the values that work best in different cases for different variables. To investigate the parameter spread around the optimal values, two types of statistics were calculated: 1. Variable-dependent statistics: Mean Error (ME), Root-Mean-Square error (RMS), and Spearman Rank-Correlation coefficient (SRC) (e.g., Wilks, 1995). These are used to in-vestigate the quality of prediction of different variables. 2. Variable-independent combined statistics, defining common measures of goodness of fit when considering all variables together: combined Mean Error (CME) , combined Root-Mean-Square error (CRMS), and combined Spearman Rank-Correlation coefficient (CSRC). These statistics, defined below, are used during calibration in order to find the optimal set of parameters for all predicted variables, where the greatest weight will be assigned to C R M S , then to C S R C , and the lowest to C M E . Each of the predicted variables spans over different ranges of values. Hence, a linear function is used to map each forecast variable into a predefined common interval in order to create variable-independent combined statistics C M E and C R M S . First, for each L E S variable the minimum and maximum values are determined = mm([E]LEs), Lf = max([E]LES) (5.1) Then, the linear-function coefficients are defined as where a has inverse units of E and b is dimensionless. The mapping has the following form to create a scaled dimensionless variable [E]s: [E}s = as[E]+bE (5.3) 62 Next, the mapping is applied to both LES and P T M variables. This new set of scaled variables has a common range, exactly (0,1) for L E S variables. The P T M variables can extend outside of this range due to prediction error. These scaled variables are then used to calculate scaled M E and R M S for each one separately. Finally, C M E and C R M S are taken as an average of all scaled M E and R M S , respectively, over all variables. The mapping procedure is not necessary in case of SRC, as it already includes ranking of correlated vectors, and C S R C can be readily taken as an average of SRC over all variables. The mathematical and physical properties of tendency Equations (4.81) and (4.82) are utilized to simplify the calibration. First, these equations can be rewritten in a form that emphasizes the K parameters: where 'Forcing' includes all the parts of the parameterization that represent forcing by mean-flow-instabilities, 'CascadeGain' is the T K E gain from energy-cascade, 'CascadeLoss' is T K E loss to energy-cascade, and 'Decay' is the asymmetry-decay term. Next, the K parameters can be rearranged: The ratios of K parameters in these tendency equations define the relative balance between different terms and properties of steady-state solutions. These characteristics can be explored by varying the X-parameter ratios. Then, to investigate the response of the model to external forcings, one can vary KF or KF directly, with the requirement of changing the other K parameters in such a way that the K ratios are constant (and steady-state solutions are preserved). These procedures allow one to investigate 'a T-tendency response' and 'a S-tendency response' to changes in KF or KF, respectively. However, all canonical cases that are investigated in this work are in quasi-equilibrium state, and a strong tendency response to either specific-force parameter is not expected. For each studied L E S case, the optimal parameters are determined in the following way. First, five values of each investigated parameter are chosen covering the entire initial parameter range'. Then the P T M model is run repeatedly with all possible combinations of those parameters (e.g., for the case with three investigated parameters it results in 125 independent P T M runs). Then the statistics are calculated and used to determine the subrange of each parameter that gives the best results. This procedure is then repeated focusing on those subranges, until it converges towards the optimal values of the investigated parameters. 5.4 Sensi t iv i ty case-studies The behavior of the P T M for each L E S canonical case, and the optimal parameter-space values are discussed in this section. The predicted profiles of forecasted mean and turbulence variables are plotted for each case. Also, the matrices of [TiS], [AS] (for convective cases), [Y], and [C] are shown. [US] and [AS] are used in place of [T] and [S] due to better graphical interpretation. 5.4.1 EKM case - idealized Ekman layer Only six predicted variables are important for this canonical case: U, V, B, and their turbulent fluxes. The tendency equations are affected by the smallest set of P T M parameters: KF, KjG, and K$L. The best results were obtained for KF = 1.0, K^G = 20., and K^L = 17.0. dtT = Kj Forcing + KjG CascadeGain — KjL CascadeLoss b\S = KF Forcing - K£ Decay (5.4) CascadeGain — CascadeLoss (5.5) 63 \rvs\ m [C] 1 J.J 2.36 J.D 2.36 d[km] 1.21 0.07 % X 1 •° 1.21 0.07 \ 3.5 2.36 1.21 0.07 ' 3.5 2,36 1.21 0.07 ' 3.5 2.36 1.21 0.07 s[km] s|kml s[km] Figure 5.1: P T M matrices for the E K M case: [715], [Y], and [C], where s is the height of a source layer, and d is the height of a destination layer. The following shading/contouring schemes are used to display the matrices for all LES-database cases. Ten equally spaced shading levels are used for [715] and [AS] with a delimiter every 0.1 (i.e. 0.1,...,0.9) of difference between the maximum and minimum values. Eleven shading levels are used for [Y] and [C], and delimiters are set to l / 2 n , with n = 1,. . . , 10, of the maximum value (i.e., 0.5, 0.25, 0.125, 0.0625, 0.03125, 0.015625, 0.0078125, 0.00390625, 0.00195313, and 0.000976563 of the maximum value) in order to show better their properties for short integration time-steps. The shading darkness increases with the value, from no shading for small values to black for large ones. KCG The model was sensitive to the parameter ratios -T?CT (tested in the range from 0.5 to 2.0) K C L K C L and -rFp- (tested in the range from 1.0 to 20.0). The value of - f a r = 17.0 was the largest one for Kj Kj which the P T M showed consistent improvement of the statistics with increasing -far. Further Kj K C L K C G increase of -far resulted in scatter of the statistics around the optimal -jfaj; ratio but no definite improvement in the prediction was seen. Two types of consistent model response to the variation of parameters were observed. First, K C G K C L decreasing jfaj; and increasing -far had similar effects; i.e., less intense mixing in the upper part of the P B L , buildup of inertial oscillations and increased errors in mean-wind profiles, and too K C L K C G low a concentration of B. Second, increasing -far requires increased 3 ^ to produce the best Kj KJ results for a given -far . This implies that the better results are obtained when the contribution T from small-scale turbulent transport is increased relative to the contribution from large-scale transport. The T-tendency response to K- (tested in the range from 0.001 to 2.0) was not significant in the subrange of 0.2 — 2.0. The values lower than approximately 0.2 led to increased errors in the upper part of the P B L ; i.e., buildup of inertial oscillations and increased errors in the mean-wind profiles, as well as insufficient mixing and too low a concentration of B. The model was not able to reach a quasi-steady state during the simulation for values below 0.01. The matrices are shown in Figure 5.1. [715] shows maxima of approximately 0.8[m/s] at the edges of the matrix due to a strong influence of nonlocal shear-forcing associated with the bottom layer of the model. The matrix extends throughout the full depth of the model domain, due to unconstrained development of turbulence. However, it was verified that [T] converges to a quasi-steady state, probably due to the limited domain-depth. [Y] shows non-zero values throughout most of the depth of the model domain, indicating that mixing takes place also in the upper part, even though the turbulent fluxes are effectively zero there due to constant values of mixed quantities. Both [Y] and [C] indicate the strongest mixing is taking place near the ground, as observed in nature. 04 0 0.01 0.02 0.03 0 0.0002 0.0005 0.0008 B[ kg/kg | w5[kg/kg-m/s] Figure 5.2: Vertical profiles for the E K M case: P T M (solid line) vs. L E S (dashed line). 65 [TVS] [Y] s t k m l s[km] s [ k m ] Figure 5.3: P T M matrices for the OOWC case: [TVS], [Y], and [C]. See caption in Figure 5.1 for plotting details. The vertical profiles are shown in Figure 5.2. The turbulent-flux profiles have magnitudes that are larger in the bottom part of P B L , and smaller in the upper part, than L E S results. The indication of inertial oscillations are visible in the upper parts of profiles for U and V, but are effectively absent in L E S profiles. In summary, even with the optimal parameter values, P T M exhibits mixing that is relatively too strong in the lower part of P B L , compared to mixing in the upper part. 5.4.2 OOWC case - neutral boundary layer The complete set of predicted variables is important for this case, which adds dependence on critical, R C , and terminal, RT, Richardson numbers. Given the weak capping inversion for the case, the optimum Richardson-number values were RC = RT = 0.35. Even though the L E S database represents mostly cases with strong capping inversion (~ 8[-K"]), the value of RC = 0.35 was chosen for the main focus of studies here, because such a strong inversion is less widespread in the real atmosphere (e.g., occurs in southwestern California, but not usually in the rest of North America). The sensitivity to both RC (tested in the range from 0.1 to 1.0) and RT (tested in the range from RC to 1.0) was significant. Increasing the value of RC lead to increased entrainment in the upper P B L . The hysteresis effect (of turbulence onset and termination) had a negative influence on the results. Increasing the value of RT above RC led to creation of a well-mixed microlayer in the upper part of the P B L ; the larger the difference RT — R C , the deeper the microlayer. The T-tendency response to KF (tested in the range from 0.001 to 2.0) was not significant in the range of 0.2 — 2.0. Decreasing the values below 0.2 lead to shift in the peaks of the wind profiles towards the ground, and decrease of the entrainment in the capping inversion. The model was not able to reach a quasi-steady state during the simulation for values below 0.01. The matrices are shown in Figure 5.3. [TVS] shows maxima of approximately 1.3[m/s] at the edges of the matrix due to strong influence of shear-forcing relative to the bottom layer of the model. Non-zero values of [Y] are limited to the depth of the P B L . Both [Y] and [C] indicate the strongest mixing is taking place near the ground. The vertical profiles are shown in Figure5.4. Despite relatively good prediction of wu, the profiles of U are too close to linear compared with L E S profiles. The V profile does not contain a relatively well-mixed middle part as predicted by L E S . Also, the magnitude of wv is too strong both in the lower and upper parts of the P B L . 0 , B, and C are predicted relatively well, wb has a trend opposite to that predicted by L E S ; i.e., the magnitude increases with height up to the capping inversion, most probably due to excessive mixing in the lower P B L . In summary, P T M shows relatively too much mixing in the lower part of the P B L compared to the upper part. 66 Figure 5.4: Vertical profiles for the OOWC case: P T M (solid line) vs. L E S (dashed line), stands for arbitrary units. 67 [TVS] [Y] 1.38 7 3 0.71 0.05 1.38 0.71 .86 1.19 0.52 s[km] 0.05 .1 1.38 0.71 1.86 1.19 0.52 s[km] 0.05 1.86 1.19 0.52 s |km] 1.38 0.71 0.05 F i gure 5.5: P T M matrices for the 24F case: [715], L415], [Y], and [C]. See caption in Figure 5.1 for plotting details. 5.4.3 24F case - convective boundary layer Only six predicted variables are important for this sensitivity study: O, B, C, and their turbulent fluxes. This case adds parameters related to asymmetry prediction in the tendency equations: Sw, KF, and K§ . The best results were obtained for Sw = 0.0, K( = 1.0, and Kg = 0.8. With Sw = 0., the P T M does not allow for convective penetration, because the asymmetric force does not contribute to generation of T K E . As a result, the depth of the convective unstable layer is equivalent to Stull's nonlocal static-stability definition. Both VTTE-wise and parcel-wise models, presented in Section 4.2.1 of Chapter 4, are equivalent to each other for this value of Sw. The model was strongly sensitive to changes of -$r (tested in the range from 0.1 to 10.0), and K D weakly to Sw (tested in the range from 0.0 to 1.0). The sensitivity to was most noticeable in the turbulent fluxes: too strong in the upper P B L and too weak in the lower P B L for small ratio values, to too weak in the entire P B L for large values. The T-tendency response to KF (tested in the range from 0.001 to 2.0) was not significant in the range 0.5 — 1.5. Decreasing the values below 0.2 led to decreased entrainment in the capping inversion. The model was not able to reach a quasi-steady state during the simulation for values below 0.01. The S-tendency response to KF (tested in the range from 0.001 to 2.0) was not significant in the range 0.2 — 2.0. Decreasing the values below 0.2 led to decreased entrainment in the capping inversion and strong super-adiabatic profiles of 6 in the lower P B L (i.e., behavior similar to that discussed below for a symmetric-forcing-only version of the P T M ) . The model was not able to reach a quasi-steady state during the simulation for values below 0.1. The matrices are shown in Figure 5.5. TVS shows maxima of approximately 2.9[m/s] near the edges (but not so tight at the edges as for shear-driven cases) of the matrix, indicating that turbulence is driven by forcing scales slightly smaller than the depth of the P B L . [AS] is relatively constant in the entire P B L , with a maximum of approximately 4.3[m/s]. It decreases to zero near the diagonal, indicating that the small-scale motions are not skewed. It also decreases to zero towards the upper left corner, where elements are dominated by buoyancy damping. [Y] and [C] are asymmetric, with upward transport extending throughout entire P B L , and with the strongest mixing near the ground. The turbulent transport is limited to the P B L . The profiles are shown in Figure 5.6. 0 shows a super-adiabatic profile in the lower half of P B L , with relatively sharp jump in the capping inversion, compared to the L E S . This is due to insufficient mixing, which can be also seen in lack of negative entrainment flux in w9. wb and wc also show no mixing in the capping inversion, and although the magnitude of the fluxes is comparable to L E S results, the maxima are shifted towards the ground, and the depth of turbulent layer is lower than in the L E S . 08 304 1.75 1.5 1.25 N 0.75 0.5 0.25 306 0 [K] 308 310 10 B[au'| 0 0.2 0.4 0.6 0.8 1 QauJ 1.75 1.5 1.25 I 1 N 0.75 0.5 0.25 -0.05 0.05 0.1 0.15 w5[K-m/sJ 1.75 1.5 1.25 I 1 N 0.75 0.5 0.25 1.75 1.5 1.25 I ' N 0.75 0.5 0.25 -0.025 0 0.025 0.075 wb[au-m/s] 0.125 -0.01 -0.006 -0.002 0 0.002 Wc|aum/s] Figure 5.6: Vertical profiles for the 24F case: P T M (solid line) vs. L E S (dashed line), [au] stands for arbitrary units. 69 Figure 5.7: 24F case: P T M profiles for the symmetric simulation (long-dash line), and asymmetric simulation (short-dash line) vs. L E S (solid line), [au] stands for arbitrary units. Also shown are [TU?]* with symmetric forcing only, side-by-side with [TU?] including asymmetric forcing. See caption in Figure 5.1 for plotting details. 70 One may ask what is the influence of the included asymmetry on the prediction. Figure 5.7 shows the results with symmetric forcing only (K~ set to zero) versus the asymmetric P T M . Apart from insufficient mixing of B and C seen in turbulent fluxes that are too small, the largest difference is in the O profile, showing a strong super-adiabatic layer extending through lowest 75% of P B L , and with much warmer values near the ground, and much cooler values in the upper half of P B L , compared to the L E S . [TU?]* for the symmetric case is also different from [7\S] with asymmetric forcing; i.e., the strongest values develop at the edges of the matrix showing strong influence of the bottom layer on the forcing structure. Also, the values decay faster towards upper part of the P B L as a result of stronger buoyancy-damping effects there. In summary, the forcing asymmetry indeed has a strong influence on both the structure of turbulence, and on turbulent-transport. It mitigates the effect of buoyancy damping in the upper P B L , and changes the structure of [US] to being driven by both asymmetric and symmetric forcing. As the best value of Sw is equal to zero, the P T M is not capable of penetrating convection in the statistically optimal parameter-set. Summary P T M shows no significant tendency response, as was expected since all of these canonical cases represent P B L s in quasi-equilibrium. The model achieved a quasi-steady state no later than in the middle of the simulation period. It should be kept in mind that simulations of evolving stable layers could introduce more significant tendency responses. The need for a large value of KjL is interesting. It implies too strong an influence of large-scale shear forcing with decay that is relatively too small, and the large values of this parameter are necessary to contradict this effect. Nevertheless, the eddy-wise structure of [TiS] for shear-driven cases indicates that large eddies are still having a undesired strong impact on properties of turbulent transport. It is very probable that transport through the whole P B L is too strong. Adding asymmetry to convective P B L closure definitely improves the results, compared to the simulation with symmetric forcing only. However, the forcing seems to be too strong near the bottom of the P B L , and buoyancy damping is too strong in the upper P B L , resulting in a super-adiabatic O profile extending too high into P B L . No definite choice can be made between VTTE-wise and parcel-wise forcing models, since they are equivalent to each other for the optimal value of Sw. 71 Chapter 6 P T M evaluation and comparison The P T M and other turbulent P B L models are evaluated against experimental databases in this chapter. Then, the sensitivity of the P T M model to discretization is analyzed. The P T M model is run using the optimal values of its parameters that were determined in Chapter 5. 6.1 P T M comparison w i t h other models for the LES-database cases The P T M is compared with four other turbulence models for the P B L : an older implementation of T3 based on a diagnostic form of the T K E equation (called T R A N in this chapter), a nonlocal forcing closure with counter-gradient correction to the flux (called K R P O ) , a 1.5 order local closure model based on T K E prediction (called MY25), and the published L E S results. Those four models span many of the common classes of turbulence models, including those used in N W P models. A brief description of first three models is given in Appendix C, and of the L E S in Appendix B . The numerical code for first three models is attached to the L E S database ( M M M - N C A R , 2003) and readily available for comparison studies. The integration time-step, A t , and vertical resolution, Az, for the P B L models used in this section are defined based on the large-eddy turnover-time, r , and the P B L depth, Z{, given by: 1. At=±T, 2. Az — Yo . L E S case-specific r and Zi are listed in Table B.2 in Appendix B . For each case, the simulated profiles of mean variables (left column in the figures) and the corresponding turbulent fluxes (right column) are plotted. Also shown are plots of P T M matrices. Inter comparison statistics are given. The same linear mapping into a common interval, as was utilized in Chapter 5 for P T M calibration, is used here to create common comparison criteria. The C M S , C R M S , and C S R C statistics are used here to compare different models relative to the L E S . Additionally, the averages of those statistics ( A C M E , A C R M S , and A C S R C , respectively) over all cases for a given model allow overall comparison. The resulting statistics for all models and all cases, including averages, are shown in Table 6.1, which will be discussed below. 6.1.1 Neutral layer with strong capping inversion - OOSC case The modeled profiles are shown in Figure 6.1. The P T M profiles of U are too close to linear compared with the L E S profiles, despite relatively good prediction of wu. The T R A N predicts a 72 Table 6.1: Combined statistics for all compared models. Smaller values of C R M S and A C R M S , mean less difference from the L E S simulations. Higher values of C S R C and A C S R C indicate better correlation with L E S . Values of C M E and A C M E are meaningful only when C R M S and A C R M S are comparable, and then smaller absolute values indicate less difference from the L E S results. Case Statistics P T M T R A N K P R O MY25 OOSC C M E 0.00240 -0.00912 -0.01929 0.00415 C R M S 0.14310 0.13500 0.07001 0.10381 C S R C 0.78237 0.89415 0.92861 0.87973 05WC C M E -0.02114 -0.03143 0.00543 -0.00927 C R M S 0.12138 0.11240 0.12575 0.12209 C S R C 0.94233 0.95732 0.93173 0.93336 03SC C M E -0.01210 -0.02819 0.00136 -0.00092 C R M S 0.09170 0.12497 0.05429 0.07536 C S R C 0.89442 0.92297 0.93631 0.88897 05SC C M E -0.00036 -0.05587 -0.00835 -0.01427 C R M S 0.14121 0.16284 0.06505 0.07514 C S R C 0.89934 0.90162 0.93261 0.89729 24SC C M E 0.01647 -0.00262 0.02736 0.01028 C R M S 0.18102 0.16738 0.11100 0.17170 C S R C 0.80575 0.86197 0.87854 0.84839 15B C M E 0.00516 -0.00012 0.00004 -0.01245 C R M S 0.16915 0.17269 0.06882 0.10741 C S R C 0.82713 0.86945 0.90886 0.90865 24B C M E 0.00857 0.00637 0.00154 -0.00209 C R M S 0.18263 0.18666 0.08958 0.10211 C S R C 0.84296 0.87246 0.95066 0.88413 Average A C M E -0.00014 -0.01782 0.00116 -0.00351 A C R M S 0.14717 0.15171 0.08351 0.10823 A C S R C 0.85633 0.89713 0.92390 0.89150 73 Figure 6.1: Vertical profiles for the OOSC case: P T M (solid line) vs. L E S (dotted line), T R A N (single-dotted dashed line), K P R O (double-dotted dashed line), and MY25 (triple-dotted dashed line), [au] stands for arbitrary units. 74 [TVS] |Y] [Cl 0.65 0.65 v 1 1 0.33 0.02 V I • ° 0.33 0.02 0.87 0.56 0.25 0.87 0.56 0.25 0.87 0.56 0.25 s[kmj slkmj s[km] Figure 6.2: P T M matrices for the OOSC case: 715, Y , and C, where s is the height of a source layer, and d is the height of a destination layer. The following shading/contouring schemes are used to display the matrices for all LES-database cases. Ten equally spaced shading levels are used for [715] and [AS] with a delimiter every 0.1 (i.e. 0.1,...,0.9) of difference between the maximum and minimum values. Eleven shading levels are used for [Y] and [C], and delimiters are set to 1/2", with n = 1,. . . , 10, of the maximum value (i.e., 0.5, 0.25, 0.125, 0.0625, 0.03125, 0.015625, 0.0078125, 0.00390625, 0.00195313, and 0.000976563 of the maximum value) in order to show better their properties for short integration time-steps. The shading darkness increases with the value, from no shading for small values to black for large ones. relatively well-mixed layer for U in the middle of the P B L , and too low values of wu. A l l models strongly underestimate V. T R A N underestimates wv, while the other models overestimate it relative to L E S , especially in the upper part of the P B L . Contrary to the other models, the P T M shows a sharp jump in the 6, B, and C profiles. P T M turbulent fluxes for those variables are consistently underestimated, with zero flux of wc, because of the lack of entrainment in the capping inversion. The model statistics show that P T M and T R A N deviate most from L E S in terms of C R M S . K P R O shows the best agreement here with LES . In terms of C S R C , P T M shows the lowest correlation with L E S results, while the other three models are relatively close to each other. The matrices are shown in Figure 6.2. [715] shows maxima of approximately l . l [m/s] at the edges of the matrix due to the strong influence of shear-forcing relative to the bottom layer of the model. Non-zero values of [Y] are limited to the depth of the P B L . Both [Y] and [C] indicate the strongest mixing is taking place near the ground. 6.1.2 Weak capping inversion - 05WC case This is the only case with a weak capping inversion. Unlike the remaining cases, it is the only one with the advantage of using the optimal value of Rc as was determined from the OOWC case; i.e., the neutral layer with weak capping inversion. The modeled profiles are compared in Figure 6.3. A l l models, including the P T M , capture the main characteristics of the P B L , except for poor prediction of U. The most distinct feature of the P T M is the lack of a well-mixed layer in the U, B, and C profiles, relative to that predicted by both L E S and T R A N (with strong negative error), but not by K P R O and MY25. The depth of the P B L for P T M is similar to that in T R A N and MY25, but less than for the other models due to less mixing in the capping inversion. P T M has strong mixing near the ground, which is noticeable in the linear profiles of U, B, and C in the bottom part of the P B L . Overall, P T M shows behavior qualitatively closer to MY25 or K P R O , than to T R A N . The model statistics show very similar results. In terms of C R M S , T R A N deviates least from the L E S , and P T M is marginally closer to the LES than K P R O and MY25 . Having comparable C R M S , P T M shows stronger negative C M E than K P R O and MY25. In terms of C S R C , K P R O 75 Figure 6.3: Vertical profiles for the 05WC case: P T M (solid line) vs. L E S (dotted line), T R A N (single-dotted dashed line), K P R O (double-dotted dashed line), and MY25 (triple-dotted dashed line), [au] stands for arbitrary units. 76 1.73 1.17 0.61 0.05 1.73 1.17 0.61 0.05 1.73 1.17 0.61 0.05 1.73 1.17 0.61 0.05 s[km] s[km| s[km] s[km] Figure 6.4: P T M matrices for the 05WC case: TVS, AS, Y , and C . See caption in Figure 6.2 for plotting details. and MY25 has lower correlation than P T M , while T R A N has the highest. The matrices from the P T M are shown in Figure 6.4. [TVS] shows maxima of approximately 2.2[m/s] at scales slightly smaller than the depth of the P B L , and with weak motions generated at the largest scales at the corners of the matrix. [AS] is relatively constant with values of approximately 3.0[m/s], and rapidly goes to zero for small-scales near the diagonal. It also decreases rapidly to zero toward the upper left corner, where buoyancy damping dominates the forcing. [Y] and [C] show the desired asymmetry of mixing with strong upward motions extending throughout most of the P B L . 6.1.3 St rong capping inversion - 03SC, 05SC, and 24SC cases Those three cases differ mostly in the strength of the surface heat-flux, although 24SC has twice as large a z, compared to the others. The vertical profiles are shown in Figures 6.5 for 03SC, 6.6 for 05SC, and 6.7 for 24SC. A l l models agree fairly well for most mean variables except U, but differences are apparent in the turbulent-flux profiles. As in the 00WC case, the P T M prediction is qualitatively closer to MY25 or K P R O , especially when profiles of U and V and their fluxes are considered. This differs from L E S and T R A N , which both create a well-mixed layer. The quality of prediction of O by the P T M model decreases with increasing strength of surface heat-flux, with a strong super-adiabatic layer throughout entire P B L for the 24SC case. u>9 also degrades with increasing surface heat-flux, with negligible entrainment flux for the 24SC case, wb characteristics for the P T M are similar to MY25 for 03SC case, and is largely underestimated for the other two cases, wc is overestimated by the P T M for 03SC, but then strongly underestimated for the other two cases. Although non-zero flux for the 24SC case indicates that a small amount of entrainment takes place, the P T M shows relatively closer agreement with the L E S model than T R A N , which does not predict entrainment at all. Again, as in the 00WC case, it seems that P T M predicts mixing that is too strong in the lower part of the P B L compared to the upper part, and that the discrepancy increases with increasing surface heat-flux. Also, mixing in the capping inversion gets relatively weaker with increasing surface heat-flux. Looking at the statistics in Table 6.1, P T M shows increasing difference from L E S with increas-ing surface heat-flux. The other models give better results for weak surface heat-flux cases, but also degrade for 24SC. P T M has consistently large deviations from the L E S , while K P R O has the least deviation, and MY25 is closer to K P R O . In terms of C R S M , P T M gives closer agreement with L E S than T R A N for weak surface heat-flux cases, but less agreement for 24SC. In terms of C S R C , P T M is slightly closer to L E S than MY25 for 00SC and 05SC, but deviates more for 24SC. Overall, P T M ' s performance is comparable to MY25, but deviates from L E S more than T R A N and K P R O . 77 Figure 6.5: Vertical profiles for the 03SC case: P T M (solid line) vs. L E S (dotted line), T R A N (single-dotted dashed line), K P R O (double-dotted dashed line), and MY25 (triple-dotted dashed line), [au] stands for arbitrary units. 78 Figure 6.6: Vertical profiles for the 05SC case: P T M (solid line) vs. L E S (dotted line), T R A N (single-dotted dashed line), K P R O (double-dotted dashed line), and MY25 (triple-dotted dashed line), [au] stands for arbitrary units. 79 Figure 6.7: Vertical profiles for the 24SC case: P T M (solid line) vs. L E S (dotted line), T R A N (single-dotted dashed line), KPRO (double-dotted dashed line), and MY25 (triple-dotted dashed line), [au] stands for arbitrary units. 80 s[km] s[km] s[km] s[km] Figure 6.8: P T M matrices for the 03SC case: DS, AS, Y , and C. See caption in Figure 6.2 for plotting details. [TVS] [ Y ] [C | 0.63 1 5 0.33 0.02 0.63 0.63 0.63 I ^ 0.33 0.02 m I 0.33 0.02 X I • ° 0.33 0.02 0.85 0.55 0.24 s[km] 0.85 0.55 0.24 s[km] 0.85 0.55 0.24 s[km] 0.85 0.55 0.24 s[km] Figure 6.9: P T M matrices for the 05SC case: DS, AS, Y , and C. The grey background in [AS] corresponds to value of zero. See caption in Figure 6.2 for plotting details. [ T V S ] [Y] [C] 1.32 ^ 0.68 0.05 i 2 P * 005 . ^^_L* 1.32 I • ° 0.68 0 0 5 1.32 0 05 1 78 1.14 0.5 sfkm] 1.78 1.14 0.5 1 s |km| 78 1.14 0.5 sfkm] 1.78 1.14 0.5 s[km] Figure 6.10: P T M matrices for the 24SC case: DS, AS, Y , and C. The grey background in [AS] corresponds to value of zero. See caption in Figure 6.2 for plotting details. 81 The matrices are shown in Figures 6.8 for 03SC, 6.9 for 05SC, and 6.10 for 24SC. [TVS] matrices show slight differences as the surface heat-flux changes. The values of maxima increase from approximately 1.3[m/s], through 1.4[m/s], to 2.7[m/s], with increasing surface heat-flux. The shapes of the matrices also change; the elements grow larger in the upper-left part of the matrices with increasing surface heat-flux. This is a result of increasing contribution from [§] to the turbulence energy-cascade, and is related to changes in [AS]. The values of maxima of [AS] are approximately 1.9 [m/s] for the cases with weaker sur-face heat-flux, and approximately 3.8[m/s] for 24SC. The shape of the matrices changes more significantly. As surface heat-flux increases, the elements grow larger in the upper-left corner of the matrices, due to the increasing depth of the super-adiabatic layer near the bottom of the model domain, and decrease steadily towards main diagonal due to strong linearity of that super-adiabatic layer. The bright spots in [AS] appear as a consequence of local intense mixing creating well-mixed microlayers, and hence have negative skewness. These microlayers in the potential-temperature profiles are due to oscillations that are visible in [Y] and [C]. Asymmetry in [AS] matrices is present only for the 03SC case, despite the weakest surface heat-flux. This peculiar behavior, especially between the seemingly similar 03SC and 05SC cases, is a result of increased values of [TVS] without a large enough increase in [AS]. It is easy to show from Equations (4.42), (4.43), and (4.47), that the difference between upward and downward velocity for given i and j is proportional to AS^^/TVS^^. Hence in the presence of strong shear-related large-scale eddies, as discussed in Section 5.4, the added turbulence-production by buoyancy forcing is enough to destroy the transport asymmetry in the large-scale eddies. On the other hand, low values of elements of [AS] at smaller scales (closer to the main diagonal) do not allow for significant asymmetry in those scales. Both those effects lead to relatively small asymmetry in [Y] and [C]. 6.1.4 Baroclinic atmosphere — 15B and 24B cases In these two cases a variation of geostrophic forcing with height (i.e., a baroclinic atmosphere), providing extra shear forcing is added. The vertical profiles are shown in Figures 6.11 for 15B and 6.12 for 24B. None of the models match the significant decrease of U with height in the upper P B L that is simulated by the L E S . T R A N produces only a slight decrease of U. P T M offers only a slight gradient reversal near the capping inversion, but only for 15B. A l l models do react to the change in forcing environment, if one compares U profiles with the 24SC case, but their response is definitely weak compared to the L E S . For the other profiles, despite the added shear in the capping inversion, T R A N does not predict entrainment at all, while P T M still underestimates entrainment and turbulent fluxes of passive scalars. Overall, all models show behavior similar to the 24SC case. P T M statistics are marginally closer to LES than T R A N in terms of C R M S , but distinctly greater than K P R O and MY25. In terms of CSRC, P T M has the worst correlation with L E S of all cases. The correlation for MY25 decreases when surface heat-flux increases, while it increases for the other three models. This trend is different when compared to the SC cases, when all models had noticeably weaker correlation for strong surface heat-flux. The matrices are shown in Figures 6.13 for 15B and 6.14 for 24B. There are no noticeable relative differences in [TVS]. Those matrices have similar properties as for the 24SC case, with maxima of 2.6[m/s] for 15B, and 2.5[m/s] for 24SC. [AS] are slightly different, as the elements of 15B start to decay toward the diagonal at larger scales than for 24B. These matrices are also similar in character to 24SC, with maxima of 3.6[m/s] for 15B and 3.5[m/s] for 24B. Interestingly, the maxima of [TVS] and [AS] are slightly decreasing with increasing surface heat-flux, and are both lower than for 24SC. Most probably, it is a consequence of added mixing due to the increased 82 Figure 6.11: Vertical profiles for the 15B case: P T M (solid line) vs. L E S (dotted line), T R A N (single-dotted dashed line), K P R O (double-dotted dashed line), and MY25 (triple-dotted dashed line), [au] stands for arbitrary units. 83 Figure 6.12: Vertical profiles for the 24B case: P T M (solid line) vs. L E S (dotted line), T R A N (single-dotted dashed line), K P R O (double-dotted dashed line), and MY25 (triple-dotted dashed line), [au] stands for arbitrary units. 84 vrvs\ [wvs\ [Y] [C] 1.27 1.27 1.27 v I • ° 0.65 0.04 w I • ° 0.65 0.04 I • ° 0.65 0.04 1.7 1.09 0.48 1.7 1.09 0.48 1.7 1.09 0.48 1.7 1.09 0.48 s ( M s[km] s[km] s[km] Figure 6.13: P T M matrices for the 15B case: TVS, AS, Y, and C. The grey background in [AS] corresponds to value of zero. See caption in Figure 6.2 for plotting details. ITVS\ [WVS] |Y] [C] 1.27 1.27 1.27 1 ? I "° 0.65 0.04 £1 I • ° 0.65 0.04 I • ° 0.65 0.04 1.7 1.09 0.48 1.7 1.09 0.48 1.7 1.09 0.48 ' 1.7 1.09 0.48 s[km] s[km] s[km] s[km] Figure 6.14: P T M matrices for the 24B case: TVS, AS, Y, and C. The grey background in [AS] corresponds to value of zero. See caption in Figure 6.2 for plotting details. shear that is causing more linear profiles of 0, but with slightly smaller gradients. [Y] and [C] show no significant asymmetry, due to strong [TVS] compared to [AS] as discussed earlier for the SC cases. 6.1.5 Summary Overall statistical comparisons of the models are shown in the bottom part of Table 6.1. On aver-age, KPRO is closest to LES among the compared models, although this is the model specifically designed with nonlocal forcing correction for convective cases in mind. In terms of ACRMS, P T M is slightly closer to LES than TRAN. Having lower ACRMS, when one can compare A C M E , P T M shows on average an LES-deviation statistic that is an order of magnitude smaller than T R A N . P T M has A C M E that is closer to LES by over a factor of 8 compared to KPRO or MY25, while ACRMS is larger by less than a factor of 2. In terms of ACSRC, P T M agrees with LES worse than the other three models. KPRO has the best correlation of all models. Examination of the vertical profiles and matrices suggest that P T M produces too much mixing at the bottom of the PBL, and too little near the top. The mixing in the capping inversion, and therefore entrainment, is too weak. The shear forcing at large scales has too strong an effect on turbulence production, and prevents the appearance of asymmetry, except for cases with weak surface heat-flux. Qualitatively, P T M shows behavior closer to the local closures (in turbulent transport sense) than to the fully nonlocal TRAN. It is, at least partly, a consequence of the assumption of limited turbulent transport-range in P T M . This suggests that the threshold-like condition of this assumption could be improved. 8 5 d[km] [TVS] d[km] [WVS] Figure 6.15: PTM matrices for the Ebert case: [TVS] and [AS], where s is the height of a source layer, and d is the height of a destination layer. Twenty equally spaced shading levels are used with a delimiter every 0.05 (i.e., 0.05,.. .,0.95) of difference between the maximum and minimum values. The shading darkness increases with the value, from no shading for small values to black for large ones. The grey background in [AS] corresponds to value of zero. 6.2 E b e r t case Ebert et al. (1989) used LES to estimate transilient matrices for a purely convective PBL. Their experimental setup and estimation method are described in Section B.2 of Appendix B. In order to compare with the Ebert case, PTM is run with the same forcing parameters. The vertical resolution was matched to that of the Ebert LES run; i.e., uniform Az = 100[m]. The time-step of the model was set to r/10 = 109.8[s]. The air-density stratification was approximated by a standard atmosphere (Holton, 1992). The PTM was run for 6r (same as the spin-up time of the Ebert LES) and at that point the [T](6r) and [§](6T) matrices were saved. These matrices were used to calculate transilient matrices, [C](6r, te), for different time intervals te into the future, corresponding to the Ebert LES estimation intervals; i.e., te = 0.02,0.2,0.6,1.0, 2.0,4.Or. This approach assumes stationarity of PBL growth, and neglects the effects of convective-transport memory. This allows the PTM matrices for different estimation intervals to be compared with the corresponding matrices from Ebert LES. PTM [715](6T) and [AS](QT) matrices are shown (in place of [T](6r) and [S](6r)) in Fig-ure 6.15. The darker areas in both contoured matrices show maxima (2.2[m/s] for [715], and 3.3[m/s] for [-415]) for elements corresponding to the largest scales of turbulent motions, with values that gradually decay toward the main diagonal (lighter shading). The oscillatory behavior and small negative elements near the main diagonal are visible in [AS]. The profiles for 0(6r) and W9(6T) are not shown here, however they have properties similar to those already discussed. Namely, O has a slightly super-adiabatic profile in the middle of the PBL, a relatively stronger super-adiabatic profile in the surface layer, and a smooth sub-adiabatic transition in the capping inversion. w9 has a linear profile with no negative entrainment flux in the upper PBL. The depth of the PBL from the PTM is approximately 1500[m], compared to 1600[m] in Ebert LES. Hence, the PBL from PTM is not as deep as from the Ebert LES. Nevertheless, the entrainment is definitely present, as verified by both B and C tracers (with the same initial 86 <3[km] PTM d[km] E b e r t F i g u r e 6.16: P T M v s . E b e r t [ C ] T O for te = 0 .02 r , w h e r e s is t h e h e i g h t o f a s o u r c e l ayer , a n d d is t h e he igh t o f a d e s t i n a t i o n layer . T w e n t y e q u a l l y s p a c e d s h a d i n g levels a re u s e d w i t h a d e l i m i t e r e v e r y 0.05 (i .e. , 0 . 0 5 , . . . ,0.95) o f difference b e t w e e n t h e m a x i m u m a n d m i n i m u m va lues . T h e s h a d i n g d a r k n e s s increases w i t h t h e va lue , f r o m n o s h a d i n g for s m a l l va lues t o b l a c k for l a rge ones . T h e g r e y b a c k g r o u n d i n E b e r t [ C ] m c o r r e s p o n d s to v a l u e o f ze ro , a n d t h e s m a l l n e g a t i v e v a l u e s o f i n E b e r t ' s m a t r i x ( p l o t t e d as b r i g h t e r t h a n b a c k g r o u n d g r e y ) a re a r t i f a c t s o f t h e n u m e r i c a l a c c u r a c y o f t h e L E S m o d e l . c o n d i t i o n s as u s e d i n t h e L E S da t abase ) , s h o w i n g n o n - z e r o t u r b u l e n t f luxes i n t h e e n t i r e P B L . P T M a n d E b e r t t r a n s i l i e n t m a t r i c e s for different va lues o f te a r e s h o w n s i d e - b y - s i d e i n F i g -ures 6.16 t h r o u g h 6.22. A l l p l o t t e d t r a n s i l i e n t m a t r i c e s were n o r m a l i z e d u s i n g t h e m a s s o f the d e s t i n a t i o n layer , a n d s y m b o l [ C ] m is u sed for s u c h n o r m a l i z e d [C] . A s c a n b e seen, P T M success-f u l l y c a p t u r e s t h e m a i n features o f the m i x e d layer , i n c l u d i n g t h e c o n t r i b u t i o n of t h e l a rge r - s i ze n o n l o c a l edd ies t o t h e t o t a l m i x i n g . I t a l so cap tu re s t h e v a r i a t i o n o f t h e T M f r o m s h o r t t i m e - s t e p ( w h e n i t i s c lose t o t h e i d e n t i t y m a t r i x ) t o longe r t i m e i n t e r v a l s ( w h e n t h e m a t r i x h a s n e a r l y u n i f o r m va lues t h r o u g h o u t t h e t u r b u l e n t d o m a i n ) . I t a l so p u t s the m i x i n g o n l y w h e r e i t b e l o n g s ( i n t h e t u r b u l e n t P B L ) a n d n o t i n t h e free a t m o s p h e r e a b o v e . T h e s e a re k e y c h a r a c t e r i s t i c s t h a t P T M c a p t u r e s as w e l l as ea r l i e r T 3 p a r a m e t e r i z a t i o n s . O n e o f t h e i n n o v a t i o n s i n P T M ( tha t is l a c k i n g i n t h e e a r l i e r t r a n s i l i e n t m o d e l s ) i s t h e a b i l i t y for P T M t o d e s c r i b e s o m e of t h e a s y m m e t r i c a spec t s of c o n v e c t i v e m i x i n g . T h e e a r l i e r T 3 p a r a m e t e r i z a t i o n s , d i s c u s s e d i n C h a p t e r 3, b y de s ign a l w a y s p r o d u c e d s y m m e t r i c m a t r i c e s , a n d t h u s c o u l d neve r d e s c r i b e a n y a s y m m e t r y . A l s o , e a r l y a t t e m p t s t o c r e a t e a s y m m e t r i c m a t r i c e s ( S t u l l et a l . , 1989) c o u l d b e m a d e t o w o r k for o n l y one t i m e i n t e r v a l te, b u t w o u l d f a i l for o t h e r te d u e t o t r a n s i l i e n t m a t r i c e s w i t h n o n p h y s i c a l (nega t ive o r i m a g i n a r y ) e l e m e n t s . W h i l e t h e P T M is a p p a r e n t l y a s t ep i n the r igh t d i r e c t i o n ( b y i n c l u d i n g a s y m m e t r y ) , a c r i t i c a l e x a m i n a t i o n o f t h e s u b t l e de t a i l s o f the T M f r o m P T M is d o n e here . T h e r e s u l t i n g findings w i l l b e use fu l for f u tu r e p a r a m e t e r i z a t i o n efforts for a s y m m e t r i c t r a n s i l i e n t m a t r i c e s . I n o r d e r t o a i d t h e a n a l y s i s , t h e t u r b u l e n c e i n v a r i a n t a n d a s y m m e t r y i n d i c e s (p re sen ted o n page 25) were c a l c u l a t e d , a n d are l i s t e d i n T a b l e 6.2. E b e r t et a l . (1989) l i s t e d a set o f i m p o r t a n t features o f t u r b u l e n t t r a n s p o r t as r e p r e s e n t e d b y t h e i r e s t i m a t e d T M s . T h e s e are d i s c u s s e d nex t , s i d e - b y s i d e w i t h t h e c o r r e s p o n d i n g b e h a v i o r o f T M s f r o m t h e P T M . 87 d[km] PTM d[km] Ebert Figure 6.17: P T M vs. Ebert [C] m for te = 0 . 2 r . The grey background in Ebert [C] m corresponds to value of zero, and the small negative values of in Ebert's matrix (plotted as brighter than background grey) are artifacts of the numerical accuracy of the L E S model. See caption in Figure 6.16 for plotting details. d[km] PTM d[km] Ebert Figure 6.18: P T M vs. Ebert [C] m for te = 0 . 5 r . See caption in Figure 6.16 for plotting details. 88 d[km] 2 . PTM d[km] Ebert 1. 1. s[km] 1. s[km] Figure 6.21: P T M vs. Ebert [C]m for te = 2.r. See caption in Figure 6.16 for plotting details. d[km] PTM d[km] Ebert 1. v v \ 2 . 1. sfkm] 2. 1. 2. 1. s[km] Figure 6.22: P T M vs. Ebert [C]m for te = 4.r. See caption in Figure 6.16 for plotting details. 90 Table 6.2: Turbulence invariant, Im, and asymmetry indices, A\ through A 4 , for transilient matrices predicted by P T M and measured by Ebert et al. (1989). The indices with values smaller than 1 0 - 1 0 were replaced with zeros. Source Ai A2 ^ 3 A^ 0.02r P T M 0.03929 0. 0. 0. 0. Ebert 0.02078 0.00018 0.19400 0.03980 -0.02205 0.2T P T M 0.51899 0.03190 3.89862 0.15110 -2.22363 Ebert 0.41989 0.06001 1.74571 0.23335 -0.98700 0.5r P T M 1.11691 0.00423 2.18945 0.03446 -0.97466 Ebert 0.99674 0.04339 3.65256 0.11345 -1.57381 0.6r P T M 1.28791 0.00149 1.63371 0.01685 -0.66870 Ebert 1.14745 0.03612 4.07280 0.08951 -1.61903 l.Or P T M 1.50667 0.00019 0.64154 0.00520 -0.24053 Ebert 1.62532 0.02091 4.83342 0.04708 -1.42011 2.0r P T M 1.61341 0. 0.05579 0.00041 -0.01873 Ebert 2.20139 0.00783 3.42922 0.01895 -0.19729 4.0r P T M 1.62158 0. 0.00093 0. -0.00022 Ebert 2.32964 0.00209 1.70261 0.00221 +0.73982 • Mixing takes place over the whole depth of the P B L for te > l.Or as indicated by non-zero elements in the entire updraft (upper-right) part of the Ebert L E S T M s . The P T M predicts deep upward transport already for te > 0.2r, indicating upward motions that are too strong for the short te. The turbulence invariant, I m , indicates that the turbulence intensity predicted by P T M for te < 1.0T is stronger than that estimated from Ebert L E S . For larger te, Im of P T M does not grow as much as for Ebert L E S , but the increase for Ebert L E S is partly artificial due to dependence of Im on the number of [C] elements in turbulent domain. Namely, the number of turbulent layers predicted by P T M remains constant, but in the Ebert L E S the number of turbulent layers is increasing due to increasing depth of P B L . • Most of the air originating near the surface is transported elsewhere in the turbulent domain for te = 2.r as indicated by the relative minimum in the lower-right corner of estimated T M . This feature is not predicted by P T M . It is most probably due to the convective transport memory effect, which cannot be reproduced by P T M due to the numerical constraints of low tensor-order of the parameterization. • Convective overturning takes place in the estimated matrices, as indicated by local maxima away from main diagonal for te = 2.r, and a ridge of large values perpendicular to the main diagonal in the lower-left part of the matrix. P T M does not properly reproduce convective overturning. That indicates improper scale-wise structure of [T]. Most probably, small-scale forcing has too strong an influence on forcing and turbulence production, relative to large-scale forcing. • The estimated matrices are asymmetric, clearly visible for te > 0.2r as the ridge of larger elements below main diagonal, and by the larger range of elements away from main diagonal in the upper-right. P T M does reproduce some qualitative features of asymmetry. Upon inspecting the matrices, one can see that asymmetry is too strong for short rje, due to 91 excessive turbulence and skewness in P T M . Then, the matrices converge to the well-mixed limit far too quickly. The ridge forming below the main diagonal, very distinct in Ebert L E S estimations, does not appear in the P T M prediction at all. This qualitative assessment is confirmed by the following analysis of asymmetry indices: 1. A l l asymmetry indices are approximately zero for te = 0.02r, because for very short time intervals the transport is too local, and significant asymmetry is removed during the [Y] balancing. 2. The first index of absolute spectra asymmetry, A\, shows that the overall asymmetry predicted by P T M is too weak compared to Ebert L E S estimations. P T M predicts stronger asymmetry for te = 0.2r, with magnitude comparable to those estimated by Ebert L E S for te = 0.6r. Then, the index decreases much faster to zero with increasing te than that from the Ebert L E S . Again, this is an unwanted byproduct of the [Y] balancing. 3. The second index of absolute spectral asymmetry, A2 (associated with larger transport-scales), shows that the asymmetry predicted by P T M is also too weak and decays too rapidly with increasing te, compared with Ebert LES estimations. The only case where P T M develops long-range asymmetry stronger than that in Ebert L E S is for te — 0.2r, and is closest to that estimated by Ebert L E S for te = 0.6r. 4. Spectral asymmetry indices ^ 3 and A 4 indicate the sign of asymmetry in addition to its intensity. Those indices show magnitudes similar to their absolute counterparts. Nevertheless, the sign of asymmetry predicted by P T M , both overall (A3) and for larger transport-scales (A4), is the same as those estimated by Ebert L E S . (The esti-mation for te = 4.Or is an exception. However, the layout of the matrix suggests that it is rather related to the growth of the P B L depth and associated with its entrain-ment.) It indicates, that the simple asymmetry model designed into P T M , is capable of reproducing some of the qualitative features of asymmetry estimated by Ebert L E S . • Numerical diffusion in the free atmosphere is clearly visible in Ebert L E S estimations for te > l . r , as indicated by non-zero elements away from main diagonal, and was explained by Ebert et al. (1989) as an artifact of the numerical schemes. P T M is free of this effect, as there is no mechanism for turbulence generation above the capping inversion, in the absence of shear. In summary, P T M does not predict some of the subtle features of transilient matrices observed in Ebert L E S estimations. However, it does qualitatively predict the general asymmetry of the transport, despite a very simple parameterization of asymmetry of turbulent motions. The intensity of mixing, and hence of turbulence, predicted by P T M is too strong compared to that estimated by Ebert's L E S . This observation supports conclusions of Section 6.1; namely, that P T M has a tendency for excessive turbulence production, resulting in turbulent mixing that is too strong. Another potential weakness of the P T M has been revealed in this section; i.e., the tendency of the [Y] balancing to excessively remove asymmetry. 6.3 Sensitivity to discretization The sensitivity of P T M to changes in vertical resolution and integration time-step is investigated in this section. Three cases, used previously for calibration, are used here; i.e., E K M , OOWC, and 24F. Those particular cases were chosen because they allow one to isolate the effect of different forcing conditions. The results obtained for these cases in the calibration process, with time-step 92 0 0.01 0.02 0.03 0 0.0002 0.0005 0.0008 Blkg/kg] wE[kg/kg-m/s] Figure 6.23: E K M case: L E S profiles (solid line), and P T M profiles for At equal r /20 (short-dash line), T/10 (medium-dash line), and T/5 (long-dash line). At = T/10 and vertical resolution Az = Zi/10, are used here as a reference state toward which the model was initially calibrated, and hence representing the 'best' prediction. The A C R M S and A C S R C statistics, but this time averaged over the three considered cases, are used to aid the sensitivity analysis. 6.3.1 Sensitivity to integration time-step The P T M was re-run for two additional time-steps, A i = T/20 and At = T /5 . Namely, the temporal resolution of the runs was either increased or decreased by factor of 2. The predicted profiles are shown in Figures 6.23 through 6.25. The T/20 and T/10 simulations are close to each other, as is visible in the profiles and confirmed by A C R M S and A C S R C statistics. The differences for the E K M case are insignificant. A small time-step influence is noticeable for the OOWC case, where turbulent fluxes, wO, wb, and wc, are smaller for T/20, indicating a decrease in entrainment through the capping inversion. A similar, yet stronger, trend in behavior of the turbulent fluxes can be seen for the 24F case. A more significant difference is noticeable in 9 for 24F, where a strong and deep super-adiabatic layer forms for T/20. The T/10 and T/5 simulations show more significant differences. The A C S R C for those two cases are similar, but A C R M S for the T/5 run is over 60% larger than that of T/10. Significant degradation of forecast quality can be seen for all cases. The largest differences are related to 93 Figure 6.24: OOWC case: L E S profiles (solid line), and P T M profiles for A i equal T / 2 0 (short-dash line), T / 1 0 (medium-dash line), and T / 5 (long-dash line), [au] stands for arbitrary units. 94 Figure 6.25: 24F case: L E S profiles (solid line), and P T M profiles for At equal r /20 (short-dash line), r /10 (medium-dash line), and T / 5 (long-dash line), [au] stands for arbitrary units. 95 changes in turbulent fluxes due to increased entrainment through the capping inversion. On the other hand, the © profile is relatively immune to increasing time-step. Among all tested models, P T M shows the strongest sensitivity to changes in integration time-step. The predictions of MY25 and K P R O are relatively independent of time-step. T R A N shows a similar, yet weaker, trend with approximately 10% difference in A C R M S between T / 2 0 and T / 5 . Recall from calibration Section 5.4, that P T M exhibited insignificant dependence on ten-dency parameters. Thus, it is believed that the threshold-like assumption of limited turbulent transport-range (Equation (4.88)) is mostly responsible for the sensitivity discussed in this sec-tion. As integration time-step increases, the mixing becomes more nonlocal. However, there is one qualitative contradiction that might not be immediately noticeable. Namely, decreasing the time-step had little effect on the E K M or OOWC cases, where shear is mostly responsible for turbulence generation, while it had strong negative effect on the shape of the © profile, and therefore on the properties of convective forcing. The response characteristics of P T M prediction to changes in time-step correspond to features of turbulent transport. Namely, the shear-driven turbulent transport is more local, and hence there is little response to decreasing the time-step. On the other hand, convection-driven turbulent transport is more nonlocal, and hence exhibits little response of © (and buoyancy forcing) to increasing time-step. 6.3.2 Sensitivity to vertical resolution The P T M model was run using additional vertical resolutions, Az = Zj/20 and Az = Zi/b. Namely, the vertical resolution of the runs was either increased or decreased by factor of 2. The predicted profiles are shown on Figures 6.26 through 6.28. The Zj/20 and Zj/10 simulations are relatively close to each other, as is visible in the profiles and confirmed using A C R M S and A C S R C . The differences for the E K M and OOWC cases are relatively insignificant. More noticeable influence can be seen in the 24F case, where there are larger differences in wb and wc. The Zj/10 and Zi/5 simulations show significant differences. A C S R C is only approximately 3% smaller for Zj /5 , but A C R M S is approximately 260% larger. A significant degradation of overall forecast quality can be seen for all cases. Again, the largest differences are related to changes in turbulent fluxes due to increased entrainment through the capping inversion. A negative value of w9 appears for the Zi/5 simulation, with a magnitude comparable to L E S results, but lower in the P B L . P T M shows the strongest sensitivity to changes in vertical resolution, among all tested models. Neither are the other models immune to those changes. A l l of the models behave similarly to P T M ; namely, they show stronger degradation of forecast quality when vertical resolution is decreased from Zj/10 to Zj /5. K P R O and MY25 have the largest increase in A C R M S of approximately 40%, while T R A N and MY25 have the largest decrease in A C S R C of approximately 6%. There are also some subtle features of the P T M response to changes in vertical resolution. First, 0 profiles become closer to the well-mixed LES profile with decreasing resolution, and w6 profiles show evidence of entrainment, as already mentioned above. Second, for the E K M case, a larger difference between P T M and LES is evident in the turbulent fluxes in the upper P B L , and tends to decrease towards the surface. Those two features suggest that another effect is more important here (in addition to the already-mentioned limited turbulent transport range). Namely, the P T M approximation of the turbulence energy-cascade also affects the overall sensi-tivity to vertical resolution via interaction of the renormalization factor (Equation (4.76)) with the structure of the forcing. It is because the number of daughter eddies decreases when the 96 Figure 6.26: E K M case: LES profiles (solid line), and P T M profiles for Az equal Zi /20 (short-dash line), Zi/'10 (medium-dash line), and Z j /5 (long-dash line). 9 7 1.75 1.5 -0 .4 -0 .3 -0 .2 -0.1 0 wu[m 2 /s 21 N 0.75 0.5 0.25 -0 .15 -0.1 -0 .05 0 W v [ m 2 / s 2 ] 301 302 303 304 305 306 -0.004 -0.002 0 0[K] i j [ t m / s ] B[au] 1.75 1.5 0 0.01 0.02 0.03 0.04 wb[au-m/s] 0 0.2 0.4 0.6 0.8 1 -0.003 -0.002 -0.001 0 C|au] Wc[aum/s] Figure 6.27: 00WC case: L E S profiles (solid line), and P T M profiles for Az equal Z j / 2 0 (short-dash line), Z j / 1 0 (medium-dash line), and Z j / 5 (long-dash line), [au] stands for arbitrary units. 98 302 304 306 308 310 © [ K ] 1.75 1.5 1.25 I 1 N 0.75 0.5 0.25 0 1.75 1.5 1.25 1 1 N 0.75 0.5 0.25 0 \ 12 B[au] •rr-/ 0.2 0.4 0.6 C[au] 0.8 -0.05 0 0.05 0.1 0.15 0.2 w5[K-m/s] 0.2 0.4 0.6 wH[au-m/s] 1.75 1.5 „ 1-25 ^ 0.75 0.5 0.25 j \ ""*-* -— \\ \ •• ~~~~ ~— •— - — ^ \\\« -0.06 -0 .04 -0 .02 Wc[aum/s] jure 6.28: 24F case: L E S profiles (solid line), and P T M profiles for Az equal Z j / 2 0 (short-dash e), Zi/10 (medium-dash line), and z;/5 (long-dash line), [au] stands for arbitrary units. 99 vertical resolution of the model decreases, and a relatively larger fraction of energy is transferred to each daughter eddy in order to mimic turbulence-energy conservation. However, the forcing by mean-flow instabilities does not have to change proportionally. This effect has the strongest influence when buoyancy damping is considered (e.g., in the upper convective P B L or in the capping inversion above neutral P B L ) : relatively more energy is transferred to the daughter eddy when vertical resolution decreases, than is lost by negative buoyancy. This leads to relative im-provement in prediction for the profiles of O and w9 in the 4F case, and increased entrainment through the capping inversion for the OOWC and 24F cases. 6 . 3 . 3 Summary The present study did not examine the sensitivity of the L E S model itself to changes in vertical or temporal resolution. The L E S resolution was roughly At « r/1000 and Az « zi/50. Those are bounded by discretization requirements of L E S modeling. It was shown in this section that the P T M parameterizations are sensitive to both integration time-step and vertical resolution. This sensitivity is stronger than for the other tested models. One likely culprit is the threshold-like assumption that limits turbulent-transport range, and causes sensitivity to changes in integration time-step and vertical resolution. Another is the renormalization factor for conservation of turbulence kinetic energy in the turbulence energy-cascade, which causes sensitivity to changes in vertical resolution. 6.4 Conclus ions Evaluation and sensitivity studies of the P T M were performed in this chapter, comparing against a set of L E S simulations comprising either neutral or convective P B L cases. The behavior of P T M was compared to that of three other P B L models typically used in numerical applications, including the former transilient turbulence model (based on a diagnostic equation for turbulence kinetic energy). The P T M model shows performance statistics comparable with the nonlocal-closure T R A N , but worse than the other two local-closure models, relative to the L E S simulations. Qualita-tively, P T M shows behavior closer to local-closure models than to T R A N , especially noticeable in the convective P B L . When running convective P B L simulations, the P T M usually creates an excessively deep super-adiabatic layer, and lacks a distinct well-mixed layer in most of the middle P B L . The prediction of asymmetry of convective motions is the major feature distinguishing P T M from other models, apart from nonlocalness of transport and forcing. It leads to distinct im-provement of prediction for convective P B L , but it does not reproduce properly the features of transilient matrices estimated in Ebert L E S experiments. A l l previously published T3 parame-terizations ignored asymmetry altogether. The P T M exhibits artificial transient oscillations in mean and turbulent-flux profiles, due to oscillations in the transilient matrix. As a result of those oscillations, transient well-mixed microlayers can form in the mean profiles and modify the small-scale properties of turbulence, such as turbulence production, intensity, or asymmetry. The behavior of the P T M as identified in this chapter leads to the following recommendations for further research into prognostic and asymmetric parameterizations. 1. Shear forcing in large scales should be reduced to reduce the amount of mixing in the entire P B L . Such reduction should be coupled with reduction in the T K E cascade-loss parameter. 100 2. Buoyancy forcing in the lower P B L should be reduced, to reduce mixing near the surface relative to mixing in the upper P B L . 3. Buoyancy damping in the upper P B L should be reduced, to allow greater mixing and entrainment in the capping inversion. 4. The threshold-like assumption in the limited turbulence-transport range approximation should be reconsidered, as it affects sensitivity of the prediction to discretization. 5. The turbulence energy-cascade approximation should be improved, to reduce sensitivity to vertical resolution. 6. Source-oriented renormalization of [Y] should be smoothed or similarly modified, to reduce vertical-profile oscillations that create sequences of well-mixed microlayers. 7. [Y] balancing should be modified, to reduce attenuation of asymmetry. The computational efficiency of P T M was not discussed so far. The highest computational cost is associated with integrating the equations for [T] and [§]. That integration time is over an order of magnitude larger than that for [Y] balancing and calculation of [C] (Equation (4.11)) together. The [Y] balancing converged in a number of steps generally smaller than the 2nd power of the number of layers in the turbulent domain. The computational costs of [Y] balancing and [C] calculation were typically similar. It is difficult at this stage to compare the computational efficiency of P T M with other tested models, because P T M was design in different programming language and was not optimized for efficiency. However, the fact that P T M uses a higher-order matrix representation of the turbulence and uses rather complex algorithms implies that P T M will be significantly more expensive than the local closure models at the same statistical-closure-order. 101 Chapter 7 Conclusions and recommendations for future work Higher-statistical-order methods of nonlocal turbulence-closure were investigated in this disser-tation, as candidates to describe unresolved eddy mixing in and above boundary layers. This research was motivated by the following deficiencies in former Transilient Turbulence Theory (T3) parameterizations: • Asymmetry: Matrix asymmetry is the nonlocal manifestation of strong narrow updrafts and wide weak downdrafts, which in local frameworks are associated with nonzero skewness statistics. However, former T3 parameterizations imposed matrix symmetry to simplify their parameterizations. This symmetry also allowed a special numerical technique to be used to convert the Mixing-Potential (MP) matrix into a mass-conserving Transilient Matrix (TM) . • Unrealistic vertical distribution of mixing in the C B L : A n exponential decay with height and very weak mixing in the entrainment zone near the top of the C B L was observed in the former T3 parameterizations. Such a mixing structure leads to unrealistic super-adiabatic profiles extending far into upper part of the C B L . • Off-diagonal versus diagonal elements of Transilient Matrices: Elements along the main diagonal of the T M represent air that is not mixed into different destination grid-cells. Although the off-diagonal elements were reasonably parameterized such that mixing is proportional to a measure of flow instability, the diagonal elements have remained an enigma. Former T3 parameterizations made risky extrapolations from the off-diagonal elements to estimate the diagonal ones. • Limited-transport range: Infinite turbulent-transport range is a common feature of nonlocal closures. Former T3 parameterizations limited the turbulent-transport range to the turbulent domain. However, they still allow the eddies to transport air to every destination in the turbulent domain, no matter how short is the integration time-step. This implies unrealistically large turbulent velocities. So motivated to advance the understanding of nonlocal turbulence parameterization, the following new concepts and models were introduced: • Turbulent-transport mean-flow state-operators: the transilient tendency-matrix as a ten-dency operator, and transilient matrix as an integral operator. 102 • Virtual turbulent-transport eddies and virtual turbulent-transport velocities. • Convective potential-energy, shear potential-energy, and P E Richardson number. • Symmetric and asymmetric components of turbulence-forcing by mean-flow instabilities. • Two alternative models for turbulence forcing by mean-flow instabilities: VTTE-wise forc-ing and par eel-wise forcing. • A prognostic model for most terms in the nonlocal turbulence kinetic energy ( T K E ) budget. • A prognostic force-decay model for the 3rd statistical moment of the vertical velocity. To aid the evaluation of these concepts, a prognostic transilient matrix model (PTM) was cre-ated as a testbed. The P T M utilizes eight dimensionless parameters in two prognostic equations. The parameter-space was investigated using the data obtained from the numerical large-eddy sim-ulations (LES) of neutral and convective planetary boundary layers (PBL) . The concepts within the P T M testbed were evaluated and compared with three other turbulence-closure models, and with two published L E S data sets. The findings from the testbed runs are summarized here with respect to the previously listed deficiencies in former T3 parameterizations: • Asymmetry: The simple parameterization of asymmetry, implemented in the P T M test-bed, qualitatively predicts the overall asymmetry of turbulent transport, as intended. Still needing improvement are the intensity of turbulence, amount of asymmetry, and convective overturning. • Unrealistic vertical distribution of mixing in the C B L : This deficiency of former parameterizations of T M s is more evident in the P T M . Namely, the intended improvement in asymmetry was offset by an unintentional deterioration of mixing. A possible remedy is further improvements to the height- and scale-wise structure of buoyancy forcing. • Off-diagonal versus diagonal elements of Transilient Matrices: The utilization of turbulent-transport mean-flow state operators removed this deficiency. • Limited transport-range: The limited transport-range approximation, while a physically justified attribute, leads to sensitivity to changes in vertical resolution as implemented in the P T M . Based on those findings, it is recommended that the following issues are ripe for future inves-tigation: • The scale structure of nonlocal shear forcing. , • The height and scale structure of nonlocal buoyancy forcing. • Parameterizations for nonlocal turbulent transport and pressure redistribution of T K E . • Properties of asymmetry of turbulent transport in convective P B L s , and ways to enhance asymmetry in TMs . • The limited transport-range of finite-velocity turbulence, and its implementation. 103 • The turbulence energy-cascade approximation in nonlocal physical space, as opposed to traditional methods in wavenumber space. • The assumption of isotropy of the virtual turbulent-transport eddies. • Numerical methods for accelerating computations. Finally, the future work might find it useful to employ the following methodologies: • Calibration of a testbed such as an improved P T M should include cases for the stable P B L . • The verification of testbed models should include observations from the real atmosphere. • Calibration and verification should be extended to the free atmosphere and to clear-air turbulence, to fully prove applicability as a turbulence-closure. 104 References Andren, A . , A . R. Brown, J . Graf, R J . Mason, C. H . Moeng, F . T. M . Nieuwstadt, and U . Schu-mann, 1994: Large-eddy simulation of a neutrally stratified boundary-layer - a comparison of 4 computer codes. Quarterly Journal of the Royal Meteorological Society, 120(520), 1457-1484. Andren, A . and C . -H . Moeng, 1993: Single-point closures in a neutrally stratified boundary-layer. Journal of the Atmospheric Sciences, 50(20), 3366-3379. Arya, S. P., 1988: Introduction to Micrometeorology. International Geophysics Series. Academic Press, 307pp. —, 1999: Air Pollution Meteorology and Dispersion. Oxford University Press, 310pp. Ayotte, K . W. , P. P. Sullivan, A . Andren, S. C. Doney, A . A . M . Holstag, W . G . Large, J . C. McWiliams, C . -H. Moeng, M . Otte, J . J . Tribbia, and J . C. Wyngaard, 1996: A n evaluation of neutral and convective planetary boundary-layer parameterizations relative to large-eddy simulations. Boundary-Layer Meteorology, 79, 131-175. Bagliani, M . , 1997a: Chiusura non locale della turbulenza dello strato limite atmosferico. Doctoral thesis, University of Torino, Italy. In Italian, 179pp. —, 1997b: Indices of asymmetry in the context of transilient turbulence theory. University of Torino, Italy. Personal communication. Berkowicz, R., 1984: Spectral methods for atmospheric diffusion modeling. Boundary-Layer Meteorology, 30, 201-220. Berkowicz, R. and L . P. Prahm, 1979: Generalization of K theory for turbulent diffusion. Part I: Spectral turbulent diffusion concept. Journal of Applied Meteorology, 18, 266-272. Bohren, C. F . and B. A . Albrecht, 1998: Atmospheric Thermodynamics. Oxford University Press, 402pp. Caughey, S. J . , J . C. Wyngaard, and J . C. Kaimal, 1979: Turbulence in the evolving stable boundary layer. Journal of the Atmospheric Sciences, 36(1051-1052). Cenedese, A . and G. Querzoli, 1997: Lagrangian statistics and transilient matrix measurements by P T V in a convective boundary layer. Measurement Science and Technology, 8(12), 1553-1561. Cotton, W . R. and R. A . Anthes, 1989: Storm and Cloud Dynamics, volume 44 of International Geophysics Series. Academic Press. Inc., San Diego, 883pp. Cuijpers, J . W . M . and A . A . M . Holtslag, 1998: Impact of skewness and nonlocal effects on scalar and buoyancy fluxes in convective boundary layer. Journal of the Atmospheric Sciences, 55,151-162. 105 Cuxart, J . , P. Bougeault, P. Lacarrere, J . Noilhan, and M . R. Soler, 1994: A comparison be-tween transilient turbulence theory and the exchange coefficient approach. Boundary-Layer Meteorology, 67(3), 251-276. Deardorff, J . W. , 1966: The countergradient heat flux in the lower atmosphere and in the labo-ratory. Journal of the Atmospheric Sciences, 23, 503-506. —, 1972: Numerical investigation of neutral and unstable planetary boundary layers. Journal of the Atmospheric Sciences, 29(1), 91-115. Eberhard, W . L . , W . R. Moninger, and G. A . Briggs, 1988: Plume dispersion in the convective boundary layer. Part I. C O N D O R S field experiment and example measurements. Journal of Applied Meteorology, 27, 599-616. Ebert, E . E . , U . Schumann, and R. B . Stull, 1989: Nonlocal turbulent mixing in the convective boundary layer evaluated from large-eddy simulation. Journal of the Atmospheric Sciences, 46(14), 2178-2207. Fiedler, B . H . , 1984: A n integral closure model for the vertical turbulent flux of a scalar in a mixed layer. Journal of the Atmospheric Sciences, 41(4), 674-680. Fiedler, B . H . and C. -H. Moeng, 1985: A practical integral closure model for mean vertical transport of a scalar in a convective boundary layer. Journal of the Atmospheric Sciences, 42(4), 359-363. Garratt, J . R., 1992: The Atmospheric Boundary Layer. Cambridge atmospheric and space science series. Cambridge University Press, 316pp. Greenhut, G. K . and S. J . S. Khalsa, 1982: Updraft and downdraft events in the atmospheric boundary layer over the equatorial pacific ocean. Journal of the Atmospheric Sciences, 39(8), 1803-1818. Holton, J . R., 1992: An Introduction to Dynamic Meteorology. Academic Press, 511pp, 3rd edition. Holtslag, A . A . M . and B. A . Boville, 1993: Local versus nonlocal boundary-layer diffusion in a global climate model. Journal of Climate, 6(10), 1825-1842. Holtslag, A . A . M . and C. -H. Moeng, 1991: Eddy diffusivity and countergradient transport in the convective atmospheric boundary layer. Journal of the Atmospheric Sciences, 48(14), 1690-1698. Inclan, M . G. , R. Forkel, R. Dlugi, and R. B . Stull, 1996: Application of transilient turbulence theory to study interactions between the atmospheric boundary layer and forest canopies. Boundary-Layer Meteorology, 97(4), 315-344. Kaimal, J . C. and J . J . Finnigan, 1994: Atmospheric Boundary Layer Flows. Oxford University Press, 289pp. Kaimal, J . C , J . C. Wyngaard, D. A . Haugen, O. R. Cote, Y . Izumi, S. J . Caughey, and C. J . Readings, 1976: Turbulence structure in the convective boundary layer. Journal of the Atmo-spheric Sciences, 33, 2152-2169. Kolmogorov, A . N . , 1941: The local structure of turbulence in incompressible viscous fluid for very large Reynolds number. Dokl. Akad. Nauk SSSR, 30, 9-13. In Russian. 106 Larson, V . E . , 1999: The relationship between the transilient matrix and the Green's function for the advection-diffusion equation. Journal of the Atmospheric Sciences, 56(14), 2447-2453. Lenschow, D. H . and B . B . Stankov, 1986: Length scales in the convective boundary-layer. Journal of the Atmospheric Sciences, 43(12), 1198-1209. Lesieur, M . , 1997: Turbulence in Fluids. Fluid Mechanics and Its Applications. Kluwer Academic Publishers, 515pp, 3rd edition. Louis, J . F. , 1979: A parametric model of vertical eddy fluxes in the atmosphere. Boundary-Layer Meteorology, 17, 187-202. Mahrt, L . , 1976: Mixed layer moisture structure. Monthly Weather Review, 104, 1403-1418. Mellor, G . L . and T. Yamada, 1974: A hierarchy of turbulence closure models for planetary boundary layer. Journal of the Atmospheric Sciences, 31, 1791-1806. —, 1977: Corrigendum. Journal of the Atmospheric Sciences, 34(9), 14821482. —, 1982: Development of a turbulence closure-model for geophysical fluid problems. Reviews of Geophysics, 20(4), 851-875. Miles, J . W. , 1961: On the stability of the heterogeneous shear flows. Journal of Fluid Mechanics, 10, 496-508. M M M - N C A R , 2003: L E S Database and Software for P B L Model Evaluation. Personal commu-nication. Modzelewski, H . and R. B . Stull, 1997: Implementation of Transilient Turbulence in an Opera-tional Diffusion Model. Technical Report T C N 95-275, Battelle, Columbus Division, 18pp. —, 1998a: One-and-a-half order closure model of the transilient turbulence matrix. In 32nd CMOS Congress: Atmosphere-Ocean Climate Variability. Halifax. —, 1998b: A prognostic parameterization of the transilient turbulence matrix for non-local clo-sure. In European Geophysical Society XXIII General Assembly, volume 16 of Annates Geo-physicae. European Geophysical Society, Nice, France, page 608. Moeng, C . -H. , 1984: A large-eddy-simulation model for the study of planetary boundary-layer turbulence. Journal of the Atmospheric Sciences, 41(13), 2052-2062. Moeng, C . -H . and U . Schumann, 1991: Composite structure of plumes in stratus-topped boundary-layers. Journal of the Atmospheric Sciences, 48(20), 2280-2291. Moeng, C . -H . and P. P. Sullivan, 1994: A comparison of shear- and buoyancy-driven planetary boundary layer flows. Journal of the Atmospheric Sciences, 51(7), 999-1022. Moeng, C . -H . and J . C. Wyngaard, 1984: Statistics of conservative scalars in the convective boundary layer. Journal of the Atmospheric Sciences, 41(21), 3161-3169. —, 1988: Spectral-analysis of large-eddy simulations of the convective boundary-layer. Journal of the Atmospheric Sciences, 45(23), 3573-3587. —, 1989: Evaluation of turbulent transport and dissipation closures in 2nd-order modeling. Journal of the Atmospheric Sciences, 46(14), 2311-2330. 107 Moncrieff, M . W . and J . S. A . Green, 1972: The propagation and transfer properties of steady convective overturning in shear. Quarterly Journal of Royal Meteorological Society, 98, 336-352. Obukhov, A . M . , 1941a: On the distribution of energy in the spectrum of turbulence flow. Dokl. Akad. Nauk SSSR, 32(1), 22-24. In Russian. —, 1941b: Spectral energy distribution in a turbulent flow. Izv. Akad. Nauk SSSR Ser. Geogr. Geofiz., 5(4-5), 453-466. In Russian. Pope, S. B . , 2000: Turbulent Flows. Cambridge University Press, 771pp. Prahm, L . P., R. Berkowicz, and O. Christensen, 1979: Generalization of k theory for turbu-lent diffusion. Part II: Spectral diffusivity model for plume dispersion. Journal of Applied Meteorology, 18, 273-282. Press, W . H . , B . P. Flannery, S. A . Teukolsky, and W . T. Vetterling, 1986: Numerical Recipes; The Art of Scientific Computing. Cambridge University Press, 818pp. Raymond, W . H . and R. B . Stull, 1990: Application of transilient turbulence theory to mesoscale numerical weather forecast. Monthly Weather Review, 118(12), 2471-2499. Richardson, L . F . , 1922: Weather Prediction by Numerical Processes. Cambridge University Press, 236pp. Sawford, B . L . and F . M . Guest, 1987: Lagrangian stochastic analysis of flux-gradient relation-ships in a convective boundary layer. Journal of the Atmospheric Sciences, 44(8), 1152-1165. Schmidt, H . and U . Schumann, 1989: Coherent structures of the convective boundary layer derived from large-eddy simulation. Journal of Fluid Mechanics, 200, 511-562. Sorbjan, Z., 1989: Structure of Atmospheric Boundary Layer. Prentice Hall, 317pp. Stull, R. B . , 1984: Transilient turbulence theory. Part I: The concept of eddy-mixing across finite distances. Journal of the Atmospheric Sciences, 41(23), 3351-3367. —, 1986: Transilient turbulence theory. Part III: Bulk dispersion rate and numerical stability. Journal of the Atmospheric Sciences, 43(1), 50-57. —, 1987: Transilient turbulence algorithms to model mixing across finite distances. Environmen-tal Software, 2(1), 4-12. —, 1988: Pollutant dispersion and mixed-layer modeling using asymmetric transilient matrices. In 8th AMS Symposium on Turbulence and Diffusion. American Meteorological Society, San Diego, C A , pages 263-266. —, 1991a: A comparison of parameterized vs measured transilient mixing coefficients for convec-tive mixed layer. Boundary-Layer Meteorology, 55 , 67-90. —, 1991b: An Introduction to Boundary Layer Meteorology. Atmospheric Science Library. Kluwer Academic Publishers, 666pp. —, 1991c: Static stability - an update. Bulletin of the American Meteorological Society, 72(10), 1521-1529. 108 —, 1993: Review of transilient turbulence theory and nonlocal mixing. Boundary-Layer Meteo-rology, 62, 21-96. —, 1994a: Improvements to transilient turbulence theory (T3). University of British Columbia, Canada. Personal communication. —, 1994b: A review of parameterizations schemes for turbulent boundary-layer processes. In Seminar on Parameterization of Sub-grid Scale Physical Processes. E C M W F , Reading, U . K . , pages 163-174. —, 1995: Parameterization of asymmetric transilient matrices for convective P B L s . In 11th AMS Symposium on Boundary Layers and Turbulence. American Meteorological Society, Charlotte, N C , pages 57-58. —, 2000: Meteorology for Scientists and Engineers, 2nd Ed.. Brooks/Cole, Pacific Grove, Calif., 502pp. Stull, R. B . and J . Bartnicki, 2000: Transport using transilient matrices. In Fox, P. A . and R. M . Kerr, editors, Geophysical and astrophysical convection, The fluid mechanics of astrophysics and geophysics. Gordon and Breach Science Publishers, pages 363-388. Stull, R. B . and A . G. M . Driedonks, 1987: Application of the transilient turbulence parameter-ization to atmospheric boundary-layer simulation. Boundary-Layer Meteorology, 40, 209-239. Stull, R. B . , E . E . Ebert, and S. Jascourt, 1989: Convective-structure memory and non-local turbulent mixing. University of British Columbia, Canada. Personal communication. Stull, R. B . and T. Hagesava, 1984: Transilient turbulence theory. Part II: Turbulent adjustment. Journal of the Atmospheric Sciences, 41(23), 3368-3379. Stull, R. B . and E . B . Kraus, 1987: The transilient model of the upper ocean. Journal of Geophysical Research, 92(C10), 10745-10755. Tennekes, H . and J . L . Lumley, 1972: A First Course in Turbulence. M I T Press, Cambridge, Mass., 300pp. Troen, I. and L . Mahrt, 1986: A simple-model of the atmospheric boundary-layer - sensitivity to surface evaporation. Boundary-Layer Meteorology, 37(1-2), 129-148. Weil, J . C , 1990: A diagnosis of the asymmetry in top-down and bottom-up diffusion using a Lagrangian stochastic model. Journal of the Atmospheric Sciences, 47(4), 501-515. Weisman, M . L . and J . B . Klemp, 1982: The dependence of numerically simulated convective storms on vertical wind shear and buoyancy. Monthly Weather Review, 110, 504-520. Wilks, D . S., 1995: Statistical Methods in the Atmospheric Sciences: An Introduction. Academic Press, 467pp. Willis, G . E . and J . W . Deardorff, 1976: A laboratory model of diffusion into the convective planetary boundary layer. Quarterly Journal of Royal Meteorological Society, 102, 427-445. —, 1978: A laboratory study of dispersion from an elevated source in a modified convective planetary boundary layer. Atmospheric Environment, 12, 1305-1311. 109 —, 1981: A laboratory study of dispersion from a source in the middle of the convective mixed layer. Atmospheric Environment, 15, 109-117. Wolfram, S., 2003: The Mathematica Book - Fifth Edition. Wolfram Media Inc., 1488pp. Wyngaard, J . C , 1987: A physical mechanism for the asymmetry in top-down and bottom-up diffusion. Journal of the Atmospheric Sciences, 44(7), 1083-1087. Wyngaard, J . C. and R. A . Brost, 1984: Top-down and bottom-up diffusion of a scalar in the convective boundary layer. Journal of the Atmospheric Sciences, 41(1), 102-112. Wyngaard, J . C. and J . C. Weil, 1990: Transport asymmetry in skewed turbulence. Physics of Fluids, 3(1), 155-162. Yamada, T., 1983: Simulations of nocturnal drainage flows by a q2l turbulence closure-model. Journal of the Atmospheric Sciences, 40(1), 91-106. Zhang, Q. and R. B . Stull, 1992: Alternative nonlocal description of boundary-layer turbulence. Journal of the Atmospheric Sciences, 49(23), 2267-2281. Zilitinkevich, S., V . M . Gryanik, V . N . Lykossov, and D. V . Mironov, 1999: Third-order transport and nonlocal turbulence closures for convective boundary layer. Journal of the Atmospheric Sciences, 56, 3463-3477. 110 Appendix A Nonlocal closures In this appendix, two major nonlocal closures are reviewed: spectral turbulent diffusivity theory, and integral turbulence closure. A third nonlocal method, transilient turbulence theory, was reviewed in Chapter 3. All of the nonlocal turbulence closures reviewed here are first-order closures, for modeling the turbulent flux divergence term in the equation dtE = —dzw^. They differ from K-theory (dzw£ = — KdzE) by introducing a nonlocal generalization of turbulent flux in integral form where D(z, Z)is a flux distribution-function. The above equation is a convolution of two functions, and is equivalent to K-theory in the functional limit of flux distribution D(z,Z) —> —K8(z,Z) where 5(z, Z) is Dirac's 8 distribution. One should apply nonlocal turbulence closures whenever there is a potential for nonlocal turbulent transport. This will occur due to two major circumstances. First, the integration time step of the numerical model exceeds the time-scales of 'local' turbulent transport. Second, the asymmetry of either the resolution of numerical/grid representation of the flow, or asymmetry of conditions in the flow, results in inconsistent/asymmetric representation of resolved scales. This commonly happens in numerical models if their vertical resolution is significantly finer than the horizontal - finer than the possible anisotropy in flow properties. Both above conditions are obligatory if one wants to speak about nonlocal turbulent advection, although only the first would be required to consider the classical nonlocal advection by the mean motion; e.g. the mean advection beyond the so-called CFL-criterion limit. Nonlocal turbulence closures are physically limited, as in any other models. They can be used only to transport those quantities of the flow that are conserved in the Lagrangian sense, in the absence of sources and sinks of the particular variable. This implies also that the integration time-step of any nonlocal turbulence closure is a function of the degree of stationarity of the flow. Nonlocal turbulence closure may be applied only to transport those structures of the flow that are not otherwise explicitly simulated by other physical processes in the model. Application of a nonlocal closure to larger micro-scales, such as in large eddy simulation (LES) models where convective coherent structures are resolved but shear-generated eddies are not, may serve as an example of misuse of a nonlocal closure parameterizing turbulent fluxes due to convective instability. Nonlocal turbulence closures, as presented below, differ in the way they utilize the concept of generalized nonlocal turbulent flux. Due to the definition of that quantity (Equation (A.l)) all of the closures based on it assume in some sense infinite transport rate of turbulent motions, because the fluxes, or rather flux distributions, do not depend (or depend in a limited way) on (A.l) 111 the integration time-step. In mathematical sense, continuous forms of flux-distribution functions lack spatial compactness. That allows eddies to contribute to nonlocal transport across infinite distances. The discrete form of T3 is an exception, although it still assumes that eddies have the ability to contribute to nonlocal transport across the entire turbulence domain, no matter how short is the integration time step of the model or how large the domain. Spectral Turbulent Diffusivity (STD) Spectral turbulent diffusivity theory (Berkowicz, 1984; Berkowicz and Prahm, 1979; Prahm et al., 1979) relates nonlocal transport in physical space to the spectral properties of eddies. It is a first-order closure in analytical form. It is an a priori model, where the spectral structure of turbulent eddies is prescribed. The spectral turbulent diffusivity model assumes that K is a function of wave number K; namely, K — K(K). Spectrally decomposing the local tendency equation, dtE = —dzw^, one obtains the spectral turbulent diffusivity equation: b\E(K) = K(K)dZtZE(K) (A.2) Integrating the above equation over all possible wave numbers leads to the tendency equation in physical space: dtE(z, t) = dzJ H(z, Z, t)dzE(Z)dZ (A.3) where a turbulent diffusivity transfer function is defined as: H{z,Z,t) = — [ K{K,t)eiK(z-z^dK (A.4) 27T J Berkowicz and Prahm (1979) suggested, among the other flux distributions, the following form of K(K) for practical applications *(«>*) = K * (A.5) 1 + B° {i^)J where KQ is the diffusivity at the long-time integration limit, nm is the wavenumber of the largest turbulent eddies, and BQ is a dimensionless constant. Notice that the formula is related to the most energetic wavelength of turbulence-kinetic-energy spectrum. Therefore, time-step dependence of K(K, t) is located in Km, and the theory does not model the evolution of turbulence, but only its effects. Notice that STD, being in fact an integro-differential scheme, has to rely on finite differencing both in space and time. The space differencing adds additional non-physical numerical diffusion to the scheme. The time differencing introduces an additional flux-limitation assumption to avoid negative tracer concentrations, or limits the validity of the scheme to short integration time-steps. S T D , in the presented form, was developed for homogeneous turbulence, K ^ f(z), although Berkowicz (1984) indicated general features of non-homogeneous extensions. One can see that, even though the turbulent diffusivity transfer function is related to nonlocal properties of the flow, the flux divergence term has still local character, and the closure is local in time. Integral Turbulence Closure (ITC) Starting from the nonlocal generalization of turbulent flux divergence (Equation (A. l ) ) , Fiedler (1984); Fiedler and Moeng (1985) derive the integral representation of this term in the form dM = ~ f R(z, Z)(E(Z, t) - E(z, t))dZ (A.6) 112 where a turbulent exchange function R(z, Z) represents the nonlocal transfer properties of turbu-lence. Fiedler (1984) gives a few fairly simple examples of the turbulent exchange matrix Rnm, a discrete equivalent of R(z, Z). The form of the function is assumed to be universal, and it is not directly related to actual structure of the flow and existing instabilities. For example, a simple symmetric formula was proposed where h is a size of turbulence domain, and w* is the Deardorff convective velocity. Notice, that actual conditions influence the turbulent exchange matrix through the turbulence scales. Similar to STD, ITC is an a priori closure, local in time, and the turbulent flux divergence has also a local character. However, it allows for inhomogeneous turbulence and nonlocal interaction in finite range of turbulent exchange. (A.7) 113 Appendix B LES simulations B . l L E S database for P B L model evaluation A large-eddy simulation (LES) model developed by Moeng (1984) was used at N C A R to create an L E S database for a P B L model evaluation. The detailed description of the L E S model, and verification of its performance, can be found in Moeng (1984) and Moeng and Wyngaard (1988). The L E S model has been used in numerous studies of atmospheric turbulent flow, and for investigation of P B L models, e.g. see Moeng and Wyngaard (1989), Moeng and Schumann (1991), Andren and Moeng (1993), and Moeng and Sullivan (1994). B . l . l LES model A n in-depth description of L E S modeling is beyond the scope of this work. In brief, the L E S model is designed to numerically integrate the full equations of motion at fine-enough grid-spacing to directly resolve the energy-containing large-eddy turbulent structures. Unresolved, smaller, sub-grid-scale turbulence is parameterized. The dynamic flow variables carried through the simulation are three velocity components, humidity, virtual potential temperature, and sub-grid-scale T K E (used to parameterize sub-grid-scale turbulence properties). Also, a Poisson equation for pressure perturbations is solved at each time-step to ensure incompressibility. In addition to tendency equations for dynamic and thermodynamic variables, the equations for concentration of two passive tracers are also solved: 1) tracer C with zero surface-flux, injected into only the free atmosphere as an initial condition, to simulate the entrainment process at the top of P B L ; and 2) tracer B , mimicking humidity, with specified (nonzero) surface-flux. The following boundary conditions are specified by the L E S modelers to enable the integration of the model equations. Periodic boundary conditions are imposed on the lateral boundaries, consistent with the assumption of horizontal homogeneity of the flow. A combination of a radiative boundary condition for pressure and vertical velocity, a constant-gradient for virtual potential temperature, zero-gradient of horizontal velocity, and zero sub-grid-scale T K E are used at the upper boundary of the L E S model. This prevents the generation of unphysical solutions at the top of the L E S model domain. The following forcing parameters are introduced by L E S modelers to simulate the real atmo-sphere. A mean pressure-gradient acting on the flow is specified using geostrophic-wind compo-nents, and provides the means to control the strength of the shear of mean horizontal-velocity field. Surface roughness-length and surface heat-flux are other external parameters used to de-fine interaction between the atmosphere and the bottom boundary. This enables simulation of different stability conditions developing in the flow, and hence different types of P B L . The flow conditions in free atmosphere (above the P B L ) are specified by the geostrophic-wind components, 114 Table B . l : External forcing parameters for LES simulations and P B L runs: Ug and Vg are the components of geostrophic wind, z$ is surface roughness-length, and w9 is the surface sensible heat-flux. dzQ is the initial gradient of potential temperature in the free atmosphere and CI indicates if strong capping inversion was imposed in the initial conditions of L E S simulation. (Ayotte et a l , 1996). Case Ug[m/s] Vg[m/s] zo[m] w6[K™} dzG[K/m] CI E K M 10. 0. 0.10 0.00 0.000 no OOWC 15. 0. 0.16 0.00 0.006 no 24F 0. 0. 0.16 0.24 0.003 yes OOSC 15. 0. 0.16 0.00 0.003 yes 05WC 15. 0. 0.16 0.05 0.006 no 03SC 15. 0. 0.16 0.03 0.003 yes 05SC 15. 0. 0.16 0.05 0.003 yes 24SC 15. 0. 0.16 0.24 0.003 yes 15B 10. 0.-20. 0.16 0.15 0.006 yes 24B 10. 0.-20. 0.16 0.24 0.006 yes and the statically-stable gradient of virtual potential temperature specified as initial conditions. In summary, the design of the LES model, as presented above, allows it to simulate dif-ferent classes of turbulent PBL-flows using a small set of forcing parameters: geostrophic-wind components, surface roughness-length, surface heat-flux, and the gradient of virtual potential temperature. This model was used by developers at N C A R to generate a database of simulated mean-flow and turbulence output for a set of archetypical P B L s , as described next. B . l . 2 LES database at N C A R The LES database developed for P B L model evaluation (Ayotte et al., 1996) consists of a set of simulations of shear and/or buoyancy driven turbulent flows in a dry (humidity-free) P B L . The simulated cases are limited to statically unstable and neutrally stratified flows, mainly because the L E S model works best for quasi-equilibrium flows containing large turbulent eddies. The stable boundary layers, subject to negative buoyancy forcing, are typically non-equilibrium flows with the energy-containing part of the turbulence spectrum shifting persistently towards small-scale eddies and decaying in time. Hence, it is difficult to maintain a relative independence of the L E S results from sub-grid-scale parameterization in such conditions. On the other hand, the neutral and convective P B L s represent a significant subset of possible P B L flows, and offer a valuable insight for evaluation of turbulence closures. A summary of the forcing conditions used to create the L E S database used in this study is given in Table B . l . (Since only the dry atmosphere is simulated, virtual potential temperature is equal to potential temperature.) The first case, labeled E K M , corresponds to a pure forced-convection flow of the so-called true Ekman layer, with no stable free atmosphere aloft and with constant geostrophic forcing (barotropic atmosphere) throughout entire P B L . Even though conditions to fully develop such a flow are never present in the real atmosphere, they offer a possibility to isolate the development of turbulence from any influence of temperature stratification. The two cases, OOWC and OOSC, correspond to typical neutral-layer flows in a barotropic atmosphere with either weak (OOWC) or strong (OOSC) capping temperature inversions separating the P B L from the free atmosphere. Case 24F is a pure free-convection case with a relatively large surface heat-flux and strong temperature inversion at the top of the P B L . Cases 05WC, 03SC, 05SC, and 115 Table B.2: Internal scales of L E S simulations (large-eddy overturning time, r , and boundary-layer depth, used to determine vertical resolution and integration time-step of P B L runs. Also given is the length, tsim, of P B L runs. (Ayotte et al., 1996). Case r[s] Zi[m] E K M 3600 1500 19800 OOWC 1596 1065 16512 24F 515 1033 3500 OOSC 905 448 4112 05WC 905 1087 6460 03SC 624 493 4000 05SC 516 480 3500 24SC 509 1019 3500 15B 576 968 3528 24B 508 1011 3528 24SC represent typical mixed-convection conditions in barotropic atmosphere with either weak (WC) or strong (SC) capping inversion, and with turbulence generated due to combined shear and buoyancy forcing. The last two cases, 15B and 25B, introduce the effect of baroclinicity (geostrophic wind changing with height) into a mixed-convection P B L . In all cases the Coriolis coefficient is fc = 0.0001[s _ 1], which corresponds roughly to the latitude of 43.29°N. A l l L E S simulations were done using a domain of 96x96x96 cells, except for E K M which used 80x80x120 cells. Uniform grid spacing was used in all directions. The (x, y, z) domain sizes were: 2x2x3.6 [km] for case E K M ; 5x5x2 [km] for cases OOWC, 05WC, 24F, and baroclinic cases 15B and 24B; 3x3x1 [km] for cases OOSC, 03SC, and 05SC. The spatial resolution of each L E S run was chosen to satisfy L E S modeling requirements; namely that energy-containing large-eddy structures are directly resolved. The integration time-step of L E S runs varied from case to case in the range of 0.2 — 1.0[s]. The simulations were integrated for approximately 15 large-eddy turnover times r , determined using surface velocity-scales, u* or w*, and boundary-layer depth Zi. T=(Z,/W* i f ^ > 0 . [ Zi/u* if w* = 0. v ' The resulting large-eddy turnover times and boundary-layer depths are given in Table B.2. B.l .3 Creation of the LES database Evaluation of any P B L turbulence-closure model requires vertical profiles of mean variables, both to initialize the P B L model and to evaluate the results of P B L simulations. The following procedure was used by N C A R to prepare a necessary data-sets from the L E S data. The L E S model was run for couple of large-eddy turnover times (approximately 5r in most of the cases) to allow the dynamic fields to reach a well-behaved quasi-steady state. At that time, the vertical profiles of mean variables were captured from L E S , by averaging L E S fields over the horizontal layers and over time. Those profiles can be used to initialize any P B L model. The L E S model was then run for a couple more large-eddy turnover times (to approximately 10r in most of the cases) to complete the simulation. A special approach was used by N C A R to prepare the L E S output for evaluating a P B L model. Based on the observation that the reliable estimates of higher-order statistical moments 116 cannot be obtained from spatial averages at the single time-step (Ayotte et al., 1996), the com-bination of temporal and spatial averaging was used during last couple of large-eddy turnover times (approximately 6r in most of the cases). First, the spatial averages for the single time step were produced and expressed as a function of non-dimensional height z/zi in order to account for the growth of P B L . Then, the temporal averages of those were created and interpolated back to the vertical model-grid. As a result, the mean vertical profiles of mean quantities and turbulent fluxes are representative for a time in the middle of the averaging period. These profiles can then be used to evaluate P B L models. The passive tracers were initialized for all cases, except E K M , at the same time when the P B L initialization profiles were generated, using the following values of B and C (in arbitrary units): 15.0 if 2 = 0. B = { 13.5 if 0. < z < l .Olzi (B.2) 3.0 if z > l.Olzi 0.0 if z = 0. C={ 0.0 if 0. < z < l .Olzi (B.3) 1.0 if Z > 1.012; The surface flux of C is set to zero during the simulations, while the surface flux of B is calculated using a surface-similarity relationships. For the E K M case, the passive tracer C is not simulated at all, while the passive tracer B is initialized at hour 3 of the simulation with a linear profile of 0.001(1 — z/500)[kg/kg] below z = 500[m] and with constant surface flux of 0.001[% m ~ 2 s - 1 ] . In summary, the L E S database for P B L model evaluation provides two sets of data for each L E S case: 1. a set of initial conditions consisting of initial mean profiles of horizontal wind components U and V, potential temperature 0 , and passive tracers B and C. 2. a set of evaluation profiles including mean profiles of U, V, 0 , B, and C, as well as mean profiles of turbulent fluxes for momentum um and WU, heat w9, and passive tracers wb and Wt. Those data sets allow one to run a single-column P B L model using forcing conditions summarized in Table B . l , and to compare the model results against those from the L E S simulations. B.2 Ebert case The T M were estimated by Ebert et al. (1989) using an L E S simulation of the free-convectiye P B L . The L E S model of Schmidt and Schumann (1989) was used in Ebert's case. The general numerical characteristics of that model are similar to those used for the L E S database, and are not discussed here. However, this LES model was developed in Germany, and is completely independent of the N C A R L E S . Added to this German L E S was the capability to simulate an arbitrary number of passive tracers, which allows the investigation of the dispersion properties of tracers released at different vertical levels in the model. The Ebert case is a dry-atmosphere simulation. The model domain consists of 80x80x24 cells of uniform size Ax,y,z = 100[m]. External forcing is limited to a surface heat-flux of w9 = 0.06[Km/s], and the surface roughness-length ZQ = 0.16. The model was initialized with a 1350[m] deep adiabatic layer in the bottom of the domain (simulating mixed layer) having potential temperature 0 = 300 and with a constant potential temperature gradient of dz@ = 0.003[A'/m] aloft. The integration time-step of the model was approximately 4.4[s]. The model 117 was run until quasi-steady state of the turbulence field was achieved at 6r, where the large-eddy overturning time is T — 1098[s]. At that time the P B L depth was Z{ = 1600[?TT.]. After achieving this quasi-steady state (at 6r), n = 24 tracers [\|>-?=1>n] were released at 6r with the initial fields defined as: Namely, each layer in the model was initialized with a different tracer. Then, the L E S model was run for another 4r in order to estimate the T M . The transilient coefficients were collected at times te = 0.02,0.2,0.6,1.0, 2.0,4.Or (counted since the tracer initialization) by using the horizontal average of tracer concentration [tyi] as an approximation of j-th column of the T M , i.e.: Thus the variation of transilient coefficients with time is simulated by dispersion of mean tracer concentration for that time. (B.4) C M ( t e ) = * i [ l l ( t e ) (B.S) 118 Appendix C P B L models used in comparison study Here is a brief description of the other P B L models that are used in this study for comparison with P T M : Transilient Turbulence ( T R A N ) The version of T3 model used here is based on the version described in Stull and Driedonks (1987) with minor coefficient adjustment from Stull (1993). This version of the model was described thoroughly in Chapter 3. Although, a newer formulation (Inclan et al., 1996) exists, the older version was chosen for two reasons: first, it was the subject of broader studies, including more versatile validation; second, it is much less computationally expensive than more recent version; hence, if higher computational cost of P T M was to be justified, it would need to perform better than this older parameterization. Nonlocal counter-gradient flux correction (KPRO) The model used here is based on a nonlocal correction to the local K-theory, see Equation (1.10), as proposed by Deardorff (1966). In the model used here, the Deardorff correction was extended by Troen and Mahrt (1986) and Holtslag and Boville (1993) to any quantity E, with perturbation part £, according to an equation: wi(z) = -KS (dzE(z) - 7 E ) ( C l ) where is an eddy diffusivity coefficient for E. j= is counter-gradient correction due to convection, parameterized as a function of surface properties and mixed-layer depth, Z{, 7 * - " . ^ (C2) where is the Deardorff convective velocity scale, w~, is another velocity scale based on the surface friction velocity scale u* and dimensionless gradients of either momentum or temperature. (u>= for passive tracers is the same as for temperature.) A second aspect of nonlocalness is achieved by dependence of K on mixed-layer depth: KE ~ w~z{l - z/zi)2 (C.3) This approach allows for a simple 1st order formulation of turbulent fluxes, yet it allows for counter-gradient transport in a convective P B L . 1.5 order closure with T K E prediction (MY25) The schematic version of the Mellor-Ya-mada 2.5-level model (1.5 closure order) was given on page 13. The version of the model used 119 here is based on the formulation of Mellor and Yamada (1982) and Yamada (1983) using a prognostic equation for T K E with down-gradient T K E transport, and eddy-diffusivity coefficients modified to include static-stability effects based on the gradient Richardson number. The numerical code for those closures is attached to the L E S database and is readily available for comparison studies. Also, these P B L models are a subset of models studied by Ayotte et al. (1996), where the more complete description of their implementation can be found. 120 Appendix D Numerical procedures and approximations The P T M model used in this work was built entirely using the Mathematica software package, and its internal programming language. D . l P T M equations step-by-step The final set of equations used in the P T M must be solved in the order listed here: 1. For calculating [T](t) and [§](£): (a) Equations (4.49) and (4.50) (b) Equations (4.51) and (4.53) (c) Equation (4.48) (d) Equations (4.58) and (4.59) (e) Equations (4.55) and (4.56) (f) If VTTE-wise forcing: Equations (4.67) and (4.69) (g) If parcel-wise forcing: Equations (4.68) and (4.70) (h) Equation (4.77) (i) Equations (4.81) and (4.82) 2. For calculating [Y](i): (a) Equation (4.83) (b) Equations (4.84), (4.85), and (4.86) (c) Equation (4.87) (d) Equation (4.88) with (4.89) (e) Equation (4.90) using algorithm in Section 4.3.2 3. For calculating [C]{t,At): Equation (4.94) 4. For calculating [C](£,Ai): Equation (4.11) 121 D . 2 N u m e r i c a l a p p r o x i m a t i o n s The complete tendency equations for mean-flow quantities combine external forcing with turbu-lence mixing. A simple process-decomposition method is used here, and the tendency equations for each separate process are applied at each integration time-step in the following order: 1. geostrophic forcing is applied to mean-wind components, 2. initial surface-layer destabilization is applied to mean-flow quantities in the bottom grid-cell (according to the procedure explained on page 26), 3. the tendency equations for T and S are integrated, and the T M is obtained, 4. the remaining part of surface-layer destabilization is applied, if necessary, 5. the turbulent fluxes can be diagnosed at this point if desired, using Equation (4.95), 6. the turbulent mixing is performed; i.e., the T M is applied to all mean-flow quantities. This approach was found to be numerically efficient, and hence practical. Next, brief descriptions of specific numerical approximations and conditions are given. Geostrophic forcing: A n explicit analytical solution of the geostrophic tendency equations can be found if one assumes constant geostrophic wind during any one time-step A t . Such a solution is used in the model, and has the following from: [U](t + At ) = [Ug] + ([U](t) - [Ug]) cos(/ c At) + ([V](t) - [Vg]) s in( / c At) [V](t + At ) = [Vg] + ([V](t) - [Vg]) cos(/ c At) - ([U](t) - [Ug]) s in ( / c At) [ U A ) where fc is the Coriolis coefficient. Surface forcings: A midpoint method (2nd order Runge-Kutta (Press et al., 1986)) is used for surface forcings in order to ensure numerical stability of the bottom boundary conditions. Integration of [T](t) and [S](t): The tendency equations ((4.81) and (4.82)) for [T](t) and [S](t) form a set of 2N2 nonlinear equations, coupled between the different scales by the turbulence energy-cascade. Solving this set of equations is not practical from a numerical point of view. However, the integration of those equations can be simplified because of the assumption that turbulence energy cascades only from larger to smaller scales. Hence, one can solve the equations by marching from the largest scales (corners of the matrix away from main diagonal) towards the smallest scales (elements next to main diagonal). In this way smaller-scale specific-forces (Equation (4.77)) can be calculated, based on larger scales, before they are used. Then, the final tendency equations are solved using the Mathematica internal solver, implementing a combination of fixed-step or time-step-adaptive methods (Wolfram, 2003); e.g., predictor-corrector Adams method with orders 1 through 12, implicit backward differentiation formulas with orders 1 through 5, and a family of implicit and explicit Runge-Kut ta methods. In addition, four numerical effects are also considered: 1. Equation (4.81) has trivial (equal to zero) solutions if initial conditions are equal to zero. Realistic initialization of that equation poses a numerical problem, especially during a transition from neutral to convective conditions in the boundary layer, because the forcing is based on integral properties of the flow. Hence, a constant fraction of a steady-state solution of that equation has been used to initialize the T K E . A value of 0.01 was found to ensure numerical stability of integration. 122 2. In order to avoid division by a small number in Equation (4.73), the denominator of the first term was replaced with max(T[/j], 0.015). Value of 0.015 corresponds to eddies with 7l5[/j] = 0.1[m/s]; i.e., eddies with velocity scales of at least an order of magnitude smaller than those typical for asymmetric convective motions in the atmosphere. 3. In order to avoid numerical effects resulting from strong non-stationarity of [T](t) and [S] (t), the average of those terms S[i,i](*) = H s M(*- A *) + s M(*)} is used in place of T[ij](t) and §[j,j](t) in Equation (4.77). 4. In order to avoid oscillations around zero, [§](£) is set to zero whenever it changes sign. V T T E properties: Equations (D.2) were also used to mitigate numerical effects of non-station-arity. Namely, T^jj(t) and §[jj](£) were substituted in Equations (4.46) and (4.47) defining V T T E properties' Mass-conservation algorithm: ETTMAX = 1 0 - 1 0 is set in the mass-conservation algorithm (see page 54). The mass and state conservation conditions of the T M (Equations (3.8) and (3.9)) were checked at every time-step, and their error did not exceed Err MAX as well. Calculating [C](t, At) from [C](i): A n approach used here to obtain T M from T T M (Equa-tion (4.11)) is based on matrix diagonalization. One may diagonalize matrix [C], use its diagonalized form to calculate the exponential function, and then use the diagonalization matrices to reverse the diagonalization process to obtain [C]. Specifically an eigen-system approach was implemented, which gave satisfactory computational performance. The pro-cedure is simple and involves four steps: 1. calculate eigenvalues and eigenvectors of [C](£) ([Evalues](t), [Evectors](t)) = EigenSystem([C](t)) 2. apply exponentiation to the eigenvalues [Evalues](t,At) = e[Evalues](t)-At 3. create a matrix, [DiagEval](t, At) with diagonal elements equal to the new eigenvalues, with off-diagonal elements equal to zero, 4. reverse the diagonalization to obtain [C](t): [C](t, At) = [Evector s]T[DiagEval](t, At) [Evectors]'1 This result can then be used with Equation (3.1). In summary, the model used in this work was found to be numerically stable in a large variety of simulated atmospheric conditions. 123 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0052396/manifest

Comment

Related Items