Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Modification of boundary layer roll dynamics by induced tropospheric gravity waves Heuff, Darlene N. 1996

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

Item Metadata

Download

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

Full Text

MODIFICATION OF BOUNDARY LAYER ROLL DYNAMICS BY INDUCED TROPOSPHERIC GRAVITY WAVES By Darlene N . Heuff B.Sc. (Physics & Mathematics) University of New Brunswick, 1988 M.Sc. (Mathematics) Carleton University, 1990 A THESIS S U B M I T T E D IN P A R T I A L F U L F I L 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 S T U D I E S D E P A R T M E N T O F G E O G R A P H Y A T M O S P H E R I C S C I E N C E P R O G R A M M E A N D I N S T I T U T E O F A P P L I E D M A T H E M A T I C S 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 BRITISH C O L U M B I A October, 1996 © Darlene N . Heuff, 1996 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of (j7f ran hy The University of British Columbia Vancouver, Canada Date Y W 9/94. DE-6 (2/88) Abstract This study investigates two-dimensional coupled interactions between the convectively driven planetary boundary layer and overlying stable layer as a possible mechanism for the formation of large aspect ratio horizontal roll convection (i.e. roll width to height ratio greater than ~4). A two-layer spectral model of the lower atmosphere is developed where the lower layer represents the convectively driven planetary boundary layer, and the upper layer represents the remainder of the troposphere. Each model domain is investigated separately but incorporates key dynamical influences from the other. In the lower domain, a modified Rayleigh-Benard convection model (which replaces the upper bounding plane with a corrugated, rigid lid of a specified wavelength and amplitude) is used to model the effects of induced gravity waves on convectively driven horizontal roll vortices. Two analytical techniques are used: a weakly nonlinear perturbation analysis; and a fully nonlinear analysis of a truncated set of equations. The resulting systems of equations are transformed to wavenumber space using a Fourier integral technique which allows interaction between all scales of motion. Linear analytical and nonlinear numerical results indicate that rolls are induced for all nonzero values of the Grashof number in response to warping of the inversion base. Results from the weakly nonlinear perturbation analysis are shown to be misleading when compared to results obtained using the fully nonlinear method; specifically with respect to the influence of the harmonics of the imposed wavenumber on roll dynamics. Although the model is able to produce rolls which are larger than those associated with the onset of convection, results from the fully nonlinear analysis indicate that the ii dominant scale of convection (X) is limited to a range given by (0.75A-p < A, < l.5Xp), where Xp is the preferred scale associated with classical Rayleigh-Benard convection. These results are robust as the band of wavenumbers containing the dominant scale is independent of both the wavenumber and amplitude of the corrugated upper boundary over a wide range of the parameters. This suggests a maximum increase in the scale of convection of 50% and therefore the proposed mechanism (as modeled) is not able to account for observed roll aspect ratios greater than ~5. The resonant mode of the stable upper layer using a simple one-layer, linear perturbation model has been found for a single set of atmospheric conditions. The mode with maximum potential for warping (and thus influence) at the inversion is determined based on energy considerations. The wavelengths corresponding to the resonant mode and dominant mode are found to be in good agreement with numerical results of Clark et al. (1986) and observations of Kuettner et al. (1987). Results from the current study support the interpretation of multiple scaling involving large-scale cloud formations presented by Balaji et al., (1993): i.e. observed large wavelength (-25 km) cloud formations (reported by LeMone and Meitin, 1984) are the result of a combination of important gravity wave and roll scales, rather than a direct result of roll vortices with horizontal scales which are comparable to observed cloud band widths. iii Table of Contents Abstract ii Table of Contents iv List of Tables ix List of Figures xii List of Abbreviations xv Partial List of Symbols xvii 1 Introduction 1 2 Background and Model Description 6 2.1 General Form of the Equations of Motion 6 2.2 Experiment and Theory Associated with Classical Rayleigh-Benard Convection 11 2.3 Field Observations of Horizontal Roll Vortices 17 2.4 Roll Formation and Modification Mechanisms 25 2.5 Physical Model Description 34 2.5.1 Assumptions and Regime of Validity 37 2.5.2 The Convectively Driven Planetary Boundary Layer 42 2.5.3 The Stable Upper Layer 45 2.5.4 Summary of Research Objectives 46 iv 3 Modified Rayleigh-Benard Convection and the Convectively Driven Planetary Boundary Layer 48 3.1 Modified Rayleigh-Benard Convection Model 48 3.2 Governing Equations in Physical Space 53 3.2.1 Nondimensional Form of the Governing Equations 55 3.3 Governing Equations in (^,r|) Space 57 3.3.1 Mapping of the Physical Space to (£,r|) Space 57 3.3.2 Analytical Techniques 60 3.3.3 Zeroth Order Equations 62 3.3.4 First Order Equations 63 3.3.5 Truncated Equations 64 3.4 Boundary Conditions 65 3.5 Fourier Integral Analysis 66 3.5.1 Zeroth Order Equations in Wavenumber Space 67 3.5.2 First Order Equations in Wavenumber Space 69 3.5.3 Truncated Equations in Wavenumber Space 72 3.6 Analytical Solutions of the Homogeneous and Forced Equations 73 3.6.1 Zeroth Order Equations 73 3.6.2 First Order Equations with a Linearized Forcing 74 3.6.3 Linearized Form of the Truncated Equations 80 3.6.4 Discussion of the Analytical Results 85 V 3.7 Numerical Analysis and Implementation 86 3.7.1 Transformation of the Analytical Equations to Numerical Representation 86 3.7.2 Numerical Solution Procedure 88 3.8 Numerical Results 92 3.8.1 Spectral Initialization 92 3.8.2 Results from the Zeroth Order Equations 99 3.8.3 Results from the First Order Equations 103 3.8.4 Results from the Truncated Equations 110 3.8.5 Summary and Discussion of Numerical Results 131 3.9 Application of the Results to the Atmosphere 134 3.9.1 An Atmospheric Test Case 134 3.9.2 Discussion of Results 144 4 The Stable Upper Layer and Induced Gravity Wave Field 147 4.1 Description of the Physical Model 147 4.2 Atmospheric Temperature and Velocity Profiles of Clark et al, 1986 .... 149 4.3 Governing Equations 151 4.4 A One-Layer Model of the Stable Upper Layer 154 4.4.1 Resonant Mode..... 157 4.4.2 Dominant Mode: Kinetic Energy Considerations 159 4.4.3 Particular Example 165 vi 4.5 A Three-Layer Model of the Stable Upper Layer 169 4.5.1 Interface and Boundary Conditions 172 4.5.2 Particular Example 173 4.6 Summary and Discussion of Results 174 5 Summary and Conclusions 176 5.1 Summary 176 5.2 Discussion 179 5.3 Extensions of Current Work 181 References 183 Appendices A Scaling Based on Atmospheric Parameters of Surface Heat Flux and the Convective Velocity Scale w* 191 A . 1 Surface Heat Flux 191 A.2 Velocity Scale 192 A. 3 Estimate of the Convective Velocity Scale for Rayleigh-Benard Convection 192 A.4 Estimate of the Convective Velocity Scale for Atmospherically Observed H R V 193 A. 5 Nondimensionalization 193 vii B Equations Governing the Real and Imaginary Components of the Zeroth Order, First Order, and Truncated Equations 194 B. 1 Zeroth Order Equations 194 B.2 First Order Equations 195 B. 3 Truncated Equations 198 C Simplifying — (V \ + eV, )*F arising from the Truncated Equations 199 D The Unstable Conductive Base State associated with the First Order Equations 201 D. 1 Comparison with Experimenatl Conditions 203 E Comparison of Analytical and Numerical Linear Results 205 E. l First Order Equations 205 E. 2 Truncated Equations 210 F Derivation of the Interface Conditions, and the Elements of the Coefficient Matrix from the Three-Layer Model Test Case 213 F. l Derivation of the Interface Conditions 213 F.2 Elements of the Coefficient Matrix 216 viii List of Tables 2.1 Convective regimes for air based on the Rayleigh number (Ra) 12 2.2 Prandtl number for various fluids 13 2.3 Summary of observed roll characteristics. (Etling and Brown, 1993) 18 2.4 Summary of roll observations presented in the literature 20 2.4 (cont.) Summary of roll observations presented in the literature 21 2.5 Observations of Kuettner et al. (1987) 24 3.1 Class I parameter values 96 3.2 Class II parameter values 97 3.3 First order equations. Numerical test parameters 103 3.4 Complete equations. Numerical test paramters I l l 3.5 Linear (N = 1) and nonlinear (N = 3) streamfuncion amplitudes, first vertical mode. Parameter values: G = 330, e = 0.01, T f = 12, and CIII 116 3.6 Linear (N = 1) and nonlinear (N = 3) streamfuncion amplitudes, first vertical mode. Parameter values: G = 500, e = 0.01, T f = 4, and CIII 116 3.7 Linear (N = 1) streamfunction amplitudes, first vertical mode. Parameter values: G = 330, e = 0.01, T f = 12, and C l l f 122 3.8 Nonlinear (N = 3) streamfunction amplitudes, first vertical mode. Parameter values: G = 330, £ = 0.01, T f = 12, and C l l f 123 ix 3.9 Linear (N = 1) streamfunction amplitudes, first vertical mode. Parameter values: G = 500, e = 0.01, T f = 4, and Cl lg 123 3.10 Nonlinear (N = 3) streamfunction amplitudes, first vertical mode. Parameter values: G = 500, e = 0.01, T f = 4, and Cl lg 123 3.11 Zeroth order equations. Numerical test parameters 136 3.12 Truncated equations. Numerical test parameters 141 4.1 Parameter values used to calculate the resonant modes. Three-layer model 173 E. 1 Analytical and numerical components of the forcing vector f,. Parameter values: G = 1, G = 1, b = l and Cla 206 E.2 Comparison of analytical and numerical results for a = b. Parameter values: G = l , (5 = 1, b = l and Cla 206 E.3 Analytical and numerical components of the forcing vector f,. Parameter values: G = 250, (5 = 1, b = 1 and Cla 207 E.4 Comparison of analytical and numerical results for a = b . Parameter values: G = 250, a = 1, b = 1 and Cla 207 E.5 Analytical and numerical components of the forcing vectors f,, ?3, f4. Parameter values: G = 330 , a = 1, b = 1 and CIII 208 E.6: Comparison of analytical and numerical results for a = b. Parameter values: G = 330 , (5 = 1, b = 1 and CIII 208 X E.7 Comparison of analytical and numerical results for (a + b) = a 0 . Parameter values: G = 330, a = 1, b = 1 and CIII 209 E.8: Comparison of analytical and numerical results for (a - b) = a 0 . Parameter values: G = 330 , a = 1, b = 1 and CIII 209 E.9: Comparison of analytical and numerical results for a = 0, b, 2b, 3b. Parameter values: G = 1, a = 1, b = 1, e = 0.01, T f = 4 and Cla 211 E.10: Comparison of analytical and numerical results for a = 0, b, 2b, 3b. Parameter values: G = 2 5 0 , o = l , b = l , e = 0.01, T f = 8 and Cla 211 E . l 1: Comparison of analytical and numerical results for a = 0, b, 2b, 3b. Parameter values: G = 330,a = l , b = l , e = 0.01, T f = 12 and CIII 212 xi List of Figures 2.1 Physical model of the troposphere 36 2.2 The convectively driven planetary boundary layer 42 2.3 The stable upper layer 45 3.1 Modified Rayleigh-Benard convection model domain 49 3.2 Flow chart outlining the analysis procedure 52 3.3 Mapping of ( X , Z ) to (£,TI) space 57 3.4 Discretization of the wavenumber spectrum 88 3.5(a) Class I initial conditions in mapped space: Streamfunction 95 3.5(b) Class I initial conditions in mapped space: Temperature 95 3.6(a) Class II initial conditions in mapped space: Streamfunction 98 3.6(b) Class II initial conditions in mapped space: Temperature 98 3.7(a) Class III initial conditions in mapped space: Temperature 100 3.7(b) Class III initial conditions in physical space: Temperature 100 3.8 Results from the first order equations in wavenumber space. Conductive regime. 105 3.9(a) Solution to the first order equations. Streamfunction 106 3.9(b) Solution to the first order equations. Temperature 107 3.10 Results from the first order equations in wavenumber space. Convective regime. 109 3.11 Wavenumber spectra. Conductive regime 113 3.12 Wavenumber spectra obtained from the complete equations. Convective regime. 115 xii 3.13 Effect of e on the solution obtained from the truncated equations 119 3.14 Comparison of results from the complete equations based on the value of the imposed wavenumber 'b \ in physical space. Convective regime 120 3.15 'Oldest' roll wavenumber development 121 3.16 Final spectra. Truncated equations test C440 125 3.17 Linear and nonlinear wavenumber development. Truncated equations 126 3.18 Final spectra for N64 and C64 130 3.19 Mean temperature profile. Classical versus modified R-B convection 131 3.20 Initial spectra. Magnitude of the streamfunction 135 3.21 Atmospheric test case. Classical R-B convection 137 3.22 Temparture in physical space 138 3.23 Mean temperature profile. Classical R-B convection 140 3.24 Final spectra. Magnitude of the streamfunction. Modified R-B convection 143 4.1 Topography of lower surface of the stable layer 148 4.2 Velocity and temperature profiles from Clark et al. (1986) 151 4.3 Estimated profiles of potential temperature and Brunt-Vaisala frequency 151 4.4 Wavelength associated with the Scorer parameter 152 4.5 One-layer model of the stable layer 155 4.6 Resonant mode. One-layer model of the stable upper layer 158 4.7 Dominant mode. One-layer model of the stable upper layer 163 4.8 Values of r | s c a l e d = [w a (1)] / a as a function of a 164 xiii 4.9 Kinetic energgy of the resonant and dominant modes 168 4.10 Mean wind and N 2 profiles used in the three-layer model 169 D. 1 Horizontal temperature gradient associated with the warped upper plate 202 D.2 Horizontal temperature at a height'd' above the lower plate 202 D.3 Experimental environment 203 xiv List of Abbreviations A C E Artie Cyclone Experiment AHOJEX Arctic Ice Dynamics Joint Experiment A M T E X Air Mass Transformation Experiment C86 Clark etal. (1986) CIa..e class I initial conditions using parameter set a .. e CIIe..e class II initial conditions using parameter set a .. e CIII class m initial conditions finite numerical method: fifth-order Runge-Kutta with finite step G A T E GARP Atlantic Tropical Experiment H R V horizontal roll vortex JAWS Joint Airport Weather Studies K87 Kuettner et al. (1987) KonTur Convection and Turbulence Experiment LHS left hand side M A S E X Mesoscale Air-Sea Exchange M I Z E X Marginal Ice Zone Experiment M E T R O M E X Metropolitan Meteorological Experiment P B L planetary boundary layer R-B Rayleigh-Benard RHS right hand side rk4 numerical method: fourth-order Runga-Kutta with finite step size STF streamfunction Temp temperature W87 Wurtele etal. (1987) xvi Partial List of Symbols coefficient matrix upper wavenumber cutoff Airy function amplitude of the real and imaginary component of the temperature for the solution of the zeroth order equations in wavenumber space, wavenumber critical wavenumber coefficient matrix Airy function amplitude of the real and imaginary component of the streamfunction for the solution of the zero'th order equations in wavenumber space, imposed wavenumber at the inversion forcing vector gravitational acceleration vertical length scale corresponding to the mean distance between the two wavenumber satisfying G = a 2 + 7 t 2 ) 3 ] / 2 a a 2 ] Grashof number [(gH 3AT) / (v 2 T m ) ] critical Grashof number [277t4 / 8a] xvii H( v F ,0) forcing function associated with the temperature equation in physical space h(\|/a k , 0 a k ) forcing function associated with the temperature equation in wavenumber space J(.,.) Jacobian K kernel K i H related to the Bessel function of the first kind of imaginary order and arguement. k vertical mode number L j ( 1 related to the Bessel function of the first kind of imaginary order and arguement. N vertical mode cutoff N 2 Brunt-Vaisala frequency P total pressure p dynamic pressure p 0 static pressure CK*?, G) forcing function associated with the streamfunction equation in physical space q(\|/a k , 0 a k ) forcing function associated with the streamfunction equation in wavenumber space R ideal gas constant Ra Rayleigh number R d ideal gas constant for dry air xviii Ri T „ , T , T m t U u u' w w' w* (v7¥) s x , X y z, Z Greek AT 6 5Q e K Richardson number temperature of the lower and upper surface respectively mean temperature of the surfaces [(T0 + T,) / 2] dimensional time variable mean wind horizontal velocity horizontal component of the perturbation velocity vertical velocity vertical component of the perturbation velocity convective velocity scale surface heat flux dimensional and nondimensionalized horizontal coordinate in physical space solution vector containing two or more of the unknowns in wavenumber space dimensional and nondimensionalized vertical coordinate in physical space difference in temperature between the two surfaces (T 0 - T,) Delta function surface heat flux amplitude of the warping of the upper surface molecular thermal difussivity xix Ke eddy thermal difussivity Km molecular thermal conductivity A, wavelength T| vertical coordinate in mapped space p m value of the density calculated at the mean temperature a Prandtl number [v / K ] o"e eddy Prandtl number [v e / K E ] x nondimensionalized time variable \% horizontal coordinate in mapped space ^(cj, r|) streamfunction in mapped space Tl) zeroth order streamfunction in mapped space y¥ l (£, r|) first order streamfunction in mapped space \\fak streamfunction in wavenumber space 0 ( £ , r|) temperature in mapped space © o r l ) zeroth order temperature in mapped space Q, r)) first order temperature in mapped space 6 A k temperature in wavenumber space v molecular viscosity v e eddy viscosity X X Chapter 1 Introduction The current study focuses primarily on a type of fluid circulation known as the horizontal roll vortex. This particular class of geophysically occurring phenomena differs from more familiar vorticular circulation such as synoptic low pressure systems, hurricanes, tornadoes, etc., not only in scale (both physical and energetic) but also in the orientation of the principal axis of rotation. While the listed phenomena rotate about an axis aligned with the vertical, as the name suggests, horizontal roll vortices rotate about axes which lie in the horizontal plane. Of particular interest is the scale (or size) of roll formation and factors which may influence scale selection. Secondary roll circulations may be identified as counter-rotating helical vortices with along axis length scales which are orders of magnitude greater than length scales associated with the cross-roll direction (which are on the order of the depth of the boundary layer). In general, horizontal roll vortices may result from: dynamic instabilities (including pure speed shear, pure turning shear, or a combination of speed and turning shear); convective instability associated with fluid in a thermally unstable environment; or a combination of dynamic and thermal instabilities (see for example Drazin and Reid, 1981). Whether the roll axis is aligned with the mean wind (longitudinal mode), perpendicular to the mean wind (transverse mode), or at some angle depends on the dominant source of flow instability. Of particular interest is the (near) longitudinal horizontal roll vortex mode (or HRV) which has been observed to occur within the planetary boundary layer (PBL) in association with cold air outbreaks (e.g. Walter, 1980) and daytime heating over homogeneous terrain (e.g. 1 Chapter 1: Introduction 2 Kuettner et al, 1987). H R V have also been observed within the urban environment (e.g. Kropfli and Kohn, 1978) and have even been proposed as a mechanism by which forest fires may be spread (Haines, 1982). There are two principal formation mechanisms which have been proposed to explain the existence of stable roll vortices within the planetary boundary layer: shear-organized convection (e.g. Asai, 1970 a,b); and the inflection point instability associated with the Ekman velocity profile (e.g. Brown, 1970) used to model the mean flow. These two mechanisms are discussed in detail in section 2.4. Results from linear perturbation studies of these mechanisms suggests the formation of rolls with aspect ratios (wavelength to depth of the boundary layer) of approximately 3 (section 2.4). However, field observations (e.g. Asai, 1966) indicate the existence of H R V with aspect ratios ranging from 1.8 to 12.7 (section 2.3). This discrepancy between linear theory and observations has resulted in a number of suggested mechanisms (presented in section 2.4) by which larger aspect ratio modes of convection may be obtained (i.e. rolls with aspect ratios greater than approximately 4). These include: latent heat effects (e.g. Skyes et al., 1988; Laufersweiler and Shirer, 1989); anisotropic eddy viscosity (e.g. Priestly, 1962); wave-wave interaction (e.g. Chang and Shirer, 1984; Mourad and Brown, 1990); and boundary layer-gravity wave interaction (Clark et al., 1986; Hauf and Clark, 1989). We shall see in section 2.4 that no one proposed mechanism has been able to account for all observations (section 2.3) of large aspect ratio roll formation. Previous studies of H R V formation have typically focused on the boundary layer in isolation from the remainder of the lower atmosphere and have therefore implicitly assumed that Chapter 1: Introduction 3 the flow dynamics are unaffected by external (i.e. non-boundary layer) influences (e.g. Chlond, 1987; Stensrud and Shirer, 1988; Laufersweiler and Shirer, 1989). However, there have been observations (e.g. Jaeckisch, 1968) which indicate that a strong gravity wave field may exist in the stable upper layer overlying the convectively driven PBL, and it has been hypothesized (Clark et al., 1986) that the induced gravity waves may influence the scale of convective activity. The current study develops a two-layer model of the lower atmosphere which is designed to investigate coupled interactions of H R V and overlying gravity waves as a possible mechanism for the formation of large aspect ratio rolls, and has been motivated by observations of Kuettner et al. (1987) and numerical results of Clark et al. (1986)). The lower layer of the two-layer model represents the convectively driven PBL, while the upper layer represents the remainder of the troposphere which is stable. Each of the two model domains are analyzed separately but have been designed to incorporate key dynamical influences from the other. This particular mechanism is unique in that it suggest the importance of large scale non-boundary layer phenomena on boundary layer dynamics and the need to consider the dynamics of the entire troposphere when modeling boundary layer processes. It is hoped that an increased understanding of the role, development, and properties of organized large eddies in the form of H R V will provide insight into the nonlinear dynamics of the convectively driven PBL, the interactions of the P B L and overlying gravity waves, as well as problems involving secondary flows which exist within a turbulent convective environment. Chapter 1: Introduction 4 This study of coupled boundary layer-gravity wave dynamics will be presented as follows: In Chapter 2, background information is presented starting with the general form of the equations governing fluid motion. A brief description of the four main analysis techniques used to study problems of convection is followed by an overview of the associated problem of classical Rayleigh-Benard convection. Next, field observations are outlined. Proposed mechanisms to account for observed large scale H R V are presented and evaluated based on results, applicability and field observations. Finally, a description of the two-layer physical model consisting of the convectively driven planetary boundary layer and an overlying stable layer, as well as research objectives are presented. Chapter 3 outlines the theoretical development of a model of the convectively driven PBL and contains the bulk of the thesis research. Focusing on a modification of the classical Rayleigh-Benard convection model (which incorporates the effects of forcing at the inversion by the induced gravity wave field) the effects of a rigid, isothermal, warped upper lid on nonlinear flow dynamics is investigated. Analytical development of the model includes a mapping of the two-dimensional physical domain leading ultimately to a transformation of the governing equations to wavenumber space. Two analytical techniques have been applied: a weakly nonlinear perturbation analysis which expands the flow variables in powers of the amplitude of the warping (e) which is assumed to be small, and a fully nonlinear analysis which truncates the mapping to first order in e. A linearized version of the resulting set of equations is investigated analytically while the complete integrodifferential equations are solved numerically. Finally, a numerical test case using atmospherically derived values of the parameters is considered and the results are compared to observations. Chapter 1: Introduction 5 The analysis of Chapter 3 is based on the assumption that the amplitude and preferred scale of the induced gravity field is known. Chapter 4 addresses these issues by examining the behavior of the overlying stable layer based on velocity and stability profiles derived from wind and temperature profiles used in numerical simulations of Clark et al. (1986). A linear perturbation analysis is applied to both a one-layer and three-layer model of the stable upper layer in order to determine the both the resonant mode(s) of the atmosphere and the mode(s) of maximum amplitude (and therefore influence) at the inversion base. Although the importance of considering the coupled dynamics of the convectively driven boundary layer and overlying stable upper layer plays a leading role in the motivation for the current work, each of these regions has been considered in isolation from one another. Both model domains have been designed to incorporate key dynamical elements from the other in hopes that the resulting flow behavior may accurately represent the coupled system. However, it is felt that the current model design represents the maximum influence that the gravity wave field may have on boundary layer dynamics. Therefore it is anticipated that a relaxing of the imposed boundary conditions which would accompany a fully dynamic coupling of the two layers, will result in a weaker forcing than is considered under the current model framework. The final chapter relates the results presented in chapters 3 and 4, thus developing a complete picture of the anticipated lower tropospheric dynamics. As such, chapter 5 summarizes the results of the study and presents suggestions for further work. Chapter 2 Background and Model Description 2.1 General Form of the Equations of Motion The behavior of an ideal Newtonian fluid is governed by the conservation of mass (2.1), momentum (2.2) and energy equations (2.3), and an equation of state (2.4): | P + V - ( p V ) = 0 (2.1) dt — + ( V - V ) V = F - - V P + - ( V - | i V ) V (2.2) at P P = ( V - K m V ) T - P ( V - V ) - u < £ (2.3) P = pRT (2.4) where V is the velocity vector with components (u,v,w), V is the gradient, (V • V) is the divergence of the velocity, F includes any external forces per unit mass, $ is a dissipation function, P is the pressure, T is the temperature, p is the density and R is the ideal gas constant. (For dry airR = R d = 287.04 m 2 s " 2 K _ 1 . ) Also, [i is the viscosity and the molecular thermal conductivity is given by K m . The specific heat at constant pressure c p is assumed to be constant. (The adaptation of the above set equations to those governing turbulent flow which incorporates an eddy viscosity and eddy conductivity, will be discussed in sections 2.5.1 and 3.9) Although equations (2.1) through (2.4) form a closed set for the six unknowns (i.e. velocity components (u, v, w), pressure, density and temperature) they are not readily solvable in PC, — + ( V - V ) T dt 6 Chapter 2: Background and Model Description 7 general. Four main groups of analysis techniques have been used to study problems of convection: linear theory; weakly nonlinear theory; spectral, quasi-spectral and finite-difference methods; and the fully nonlinear Fourier integral technique used by Getling (1983). A brief outline of the four groups are presented in order to illustrate the advantages and/or disadvantages of each of these methods. Linear Theory A linear perturbation analysis of (2.1)-(2.4) involves expanding the flow variables (f) about a given base state (f0) as f (t, x, y, z) = f0 (x, y, z) + ef, (t, x, y, z) (2.5) where f e [u,v,w,p,T,p] and the term involving f, represents the perturbation. The parameter e is assumed to be small, physically meaningful, and may represent (for example) the amplitude of the secondary circulation. Substituting (2.5) into (2.1)-(2.4) the terms of the equations are expanded and those of order e 2 are neglected. Equating coefficients of 8 to zero, the equations are separated into those governing the base state and those governing the perturbation. Using a normal mode analysis which assumes that the variables f, may be written in the form f,(t,x,y,z,) = e s t f,(x,y,z), the perturbation with the maximum growth rate (s) is determined. Application of linear theory typically assumes that this mode of maximum growth rate will initially dominate flow development. The assumption that terms involving the product of the perturbation quantities are negligible (i.e. second order terms), confines the regime of applicability of a linear perturbation Chapter 2: Background and Model Description 8 analysis to the initial stages of flow development for which the amplitude of the perturbations are infinitesimal. Weakly nonlinear theory An extension of the perturbation analysis to include the nonlinear terms is done via a weakly nonlinear approach which assumes that the flow variables may be written in the form f (t, x, y, z) = f0 (x, y, z) + efj (t, x, y, z) + e 2 f 2 (t, x, y, z)+... (2.6) As before f £ [u,v,w,p,T,p] and the parameter 8 is assumed small. Substituting (2.6) into (2.1)-(2.4) the coefficient of powers of e are equated to zero. A weakly nonlinear approach has definite advantages over a linear analysis as it retains more of the physics of (2.1)-(2.4) (i.e. terms to second order in £). However as we shall see in section 3.3.4, a weakly nonlinear analysis results in equations which are linear in respective powers of e, i.e. first order equations (denoted with a subscript '1') are linear in first order variables, second order equations ('2') are linear in second order variables etc. Thus, a weakly nonlinear approach does not allow for the exchange of energy between modes of the first (or higher) order equations, nor does it allow for the modification of the base state (f0) by the perturbation. The requirement that a physically meaningful parameter 8 exists and that this parameter is small, limits the applicability of a weakly nonlinear analysis to early stages of flow development when the amplitude of the perturbation is small. Chapter 2: Background and Model Description 9 Spectral, Quasi-spectral, and Finite-Difference Methods Spectral methods assume that both the horizontal and vertical dependence of the flow variables may be expanded using Fourier series. For example, a two-dimensional Fourier series representation of the variable f may be of the form f(t,x,z) = ]T ]T f k n(t)sin(kax)sin(n7tz) (2.7) n = l k = l where k is a horizontal mode number, a is a given wavenumber (assumed known) and n is a vertical mode number. The governing equations are then transformed to wavenumber space. The Fourier coefficients are usually solved numerically as an initial value problem. In practice, the series must be truncated at some finite value thus limiting the number of modes which are free to interact and exchange energy in a nonlinear model. Previous studies of convection have typically included only one or two vertical modes and two or three horizontal modes (e.g. Chang and Shirer, 1984). The term quasi-spectral is used to refer to a partial transformation of the flow variables with (say) the horizontal flow dependence assumed to be represented by a Fourier series while the vertical dependence is solved for explicitly in physical space, using (for example) finite-difference methods. Quasi-spectral models naturally require truncation at some finite value of the horizontal mode number. Both spectral and quasi-spectral methods typically assume that the horizontal form of the Fourier series is based on a predetermined wavenumber 'a' and its harmonics, and therefore does not allow full interaction between all scales of motion. Equations (2.1)-(2.4) may also be solved directly using numerical techniques such as finite-difference methods which approximate the governing equations with a discretized form over Chapter 2: Background and Model Description 10 a regular domain grid. Use of a finite horizontal domain requires specification of horizontal boundary conditions and periodicity over a given length scale is typically assumed. With respect to investigations of the preferred scale of convection (A, p), the main disadvantages of direct numerical simulations based on (2.1)-(2.4) is the horizontal and vertical resolution required in order to distinguish between (say) X =4 and A, =4.2 . Also, the prescription of periodicity at the horizontal boundaries requires that an integral number of rolls are maintained within the domain; restricting natural roll development. Getling (19831 The analysis technique used in the current study follows the fully nonlinear spectral methods of Getling (1983). Since the Fourier integral technique he used to study R-B convection is more general than the spectral methods discussed above, they are presented here separately. Further justification for the use of a Fourier integral method versus other methods will be discussed in section 2.5 following the presentation of the physical model under consideration, overall research objectives and results from previous studies. Getling applies a Reynolds decomposition which assumes that the flow variables may be separated into a mean (f) and a perturbation component (f ) of the form f(t, x, y, z) = f(x, y, z) + f (t, x, y, z) where as before f e [u,v,w,p,T,p]. Substituting the above into (2.1)-(2.4), the equations are separated into those governing the mean flow and those governing the temporal development of the perturbation. The nonlinear terms are retained. Since the expansion of the flow variables does not involve a small parameter such as e (cf. (2.5) and (2.6)), this method is fully nonlinear. Chapter 2: Background and Model Description 11 The horizontal flow dependence is then assumed to be represented by a Fourier integral rather than a Fourier series and the equations governing the perturbations are transformed to wavenumber space. For example, in two dimensions the Fourier integral representation of the flow variables may be of the form (c.f. 2.7). Such a representation of the flow variables, combined with a nonlinear set of governing equations, allows full nonlinear interaction between all scales as opposed to interactions which are restricted to those between a fundamental scale and its harmonics. We shall see in section 3.8.5 that this freedom of interaction between all scales accounts for the principal difference between results from the current and previous studies. 2.2 Experiment and Theory Associated with Classical Rayleigh-Benard Convection The problem of convection between two horizontal planes (i.e. classical R-B convection) is of particular interest as a modified R-B convection model is used to investigate the response of the convectively driven PBL to forcing at the inversion. Therefore, a brief overview of classical R-B convection is presented with particular focus on experimental observations and theoretical developments as they relate to the question of the preferred scale of convection. When a fluid at rest is confined between two flat plates of infinite horizontal extent, the stability of the flow is determined by thermal conditions imposed at the boundaries. If for example, the upper surface is warmer than the lower then the conductive state is stable and the fluid will remain at rest. However, if the lower surface is heated (relative to the upper surface) (2.8) n = Chapter 2: Background and Model Description 12 then there exists a finite range of temperature differences for which the conductive solution is stable. The bifurcation point may be expressed in relation to the nondimensional Rayleigh number which represents the balance between the destabilizing buoyancy force and the stabilizing processes of thermal diffusion and viscous drag. The Rayleigh number (Ra) may be expressed as a function of the separation distance between the plates (H), the vertical temperature gradient (dT / dz), viscosity (v) and thermal diffusivity (K) of the fluid, as Ra = g(dT / dz)H 4 / ( V K ) , where 'g' is the acceleration due to gravity. Above this critical value, buoyancy effects result in an induced circulation. The approximate values of Ra associated with various convection regimes (for air) ranging from time independent and well ordered (or laminar) to fully chaotic (or turbulent) are outlined in table 2.1 (Willis et al. 1972; Lipps 1976). Also, Krishnamurti (1970b, figure 9) presents a regime diagram outlining transitions between various convective states as a function of the Prandtl number a (where a = v / K represents the ratio of the kinematic viscosity to thermal diffusivity) ranging from 0.71 to infinity, based on experimental observations. The value of the Prandtl number for various fluids is presented in table 2.2. Range of Ra Convective regime 1708 < Ra <~ 5800 Ra ~ 5800 5800 < Ra < 8500 8500 < Ra 30,000 < Ra 2-D steady-state onset of small amplitude 3-D disturbances small amplitude 3-D disturbances (time-periodic) large amplitude 3-D disturbances (aperiodic) fully turbulent Table 2.1: Convective regimes for air based on the Rayleigh number (Ra). (Combined from Willis et al. 1972, and Lipps 1976). Chapter 2: Background and Model Description 13 Fluid Prandtl Number mercury water silicon air 0.025 0.71 6.8 450 Table 2.2: Prandtl number for various fluids. Linear Theory In order to investigate solutions of the classical R-B convection problem, the base state about which the perturbation analysis is applied (section 2.1) is the conductive (or null) state of the fluid. Linear perturbation analysis has been successful in predicting the critical Rayleigh number (Ra c) associated with the onset of the instability of the conductive solution as well its preferred wavelength (A.c) and associated wavenumber (a c =2n/Xc). Notably, Rayleigh (1916) solved the linear problem for the case of free-slip boundary conditions (Ra c =657.5, a c =2.21, Xc =2.83); Jeffreys (1926) investigated no-slip boundary conditions (Ra c =1708, a c=3.12, X, c=2.01); and Pellew and Southwell (1940), mixed boundary conditions. For supercritical values of the Rayleigh number (Ra) the concept that the aspect ratio of convection may be determined by maximizing the amount of heat transported by the convective circulation results in a preferred scale which decreases with increasing Ra (Malkus, 1954). Thus, linear theory predicts a decrease in wavelength as Ra increases which does not agree with experimental observations. Experimental investigations by (for example) Koschmieder Chapter 2: Background and Model Description 14 (1969) and Krishnamurti (1970a) indicate an increase in the preferred scale with an increase in Ra above critical, i.e. dA, / d(Ra) > 0 (where X is the observed wavelength of convection); increasing significantly for low Prandtl number fluids (e.g. Willis et al., 1972). Specifically, experimental results of Willis et al. (1972) indicate an increase in wavelength of -88% at R / R c =6 and -110% for R / R c = 18 for air with a = 0.71. Nonlinear Theory Schliiter et al. (1965) used a weakly nonlinear perturbation analysis (section 2.1) to investigate convection near the critical Rayleigh number. Their results indicate that only two-dimensional roll convection is stable within a subset of the linearly predicted instability range. It should be noted that the wavelengths associated with this stable subset are smaller than the critical value (i.e. X<XC). Nonlinear perturbation results of Busse (1967) at high Rayleigh numbers (with o" —> °°) indicate a stable subset of rolls in the range Ra / R a c = [1.0,13.2]. The assumption that the Prandtl number tends to infinity, removes the nonlinear term from the momentum equation (2.2), while the nonlinear term associated with the energy equation (2.3) is retained. It was found that under this condition, the subset of stable rolls contains, those with wavelengths both (slightly) larger and smaller than X, c . These theoretical results were later verified experimentally by Krishnamurti (1970 a,b) who found stable rolls within the range of Rayleigh numbers predicted by Busse (1967). Chapter 2: Background and Model Description 15 Numerical Results Numerical modeling efforts have lead to conflicting results with respect to the ability of the nonlinear exchange of energy, via the advection terms of (2.2) and (2.3), to account for observations which indicate an increase in the wavelength of convection above the critical value of the Rayleigh number (Lipps and Somerville, 1971; Getling, 1983; Chang and Shirer, 1984; Rothermal and Agee, 1986; Sykes and Henn, 1988). Lipps and Somerville (1971) considered both two- and three-dimensional convection using a finite-difference model, for Prandtl numbers of 0.71 and 6.8. They found that the two-dimensional version of their model resulted in rolls with wavelengths "k which were less than X,c while three-dimensional model results found k > Xc. This discrepancy lead them to suggest that the experimentally observed increase in wavelength is a result of three-dimensional transient processes, even if the flow is tending to a two-dimensional steady state. As discussed in section 2.1, Getling (1983) used a fully nonlinear analysis and Fourier integral techniques to investigate flow behavior over a continuous wavenumber spectrum. Getling found an increase in the preferred scale of convection as Ra increases for a = 0.1. However for values of the Prandtl number of 1 and 10 no such shift to larger wavelengths was observed and in fact a small decrease in the preferred scale was obtained. Tendency towards a single, unique preferred scale was found in contrast to the range of stable values predicted by e.g. Busse (1967). The apparent range of stable rolls observed in laboratory experiments (Krishnamurti, 1970 a,b) is attributed by Getling to the stabilizing effect of existing rolls, artificial initial roll inducement, and possible boundary effects. Also, Getling attributed the two-dimensional results of Lipps and Somerville (1971) to the initial dominance of linear processes for small amplitude disturbances Chapter 2: Background and Model Description 16 (i.e. the tendency towards smaller scales) and the stability associated with the numerical model once the flow has developed, thus inhibiting roll adjustment. Chang and Shirer (1984) used a weakly nonlinear analysis combined with a two-dimensional, seven-coefficient truncated spectral representation of (2.1)-(2.4) to study convection for small Prandtl number flows. In agreement with results of Getling (1983), they found that the preferred wavelength increased with increasing Ra for a = 0.1. Rothermal and Agee (1986) used a large aspect ratio computational domain to investigate atmospheric flow behavior using a using a quasi-spectral, two-dimensional model (section 2.1) of classical R-B convection. They found that for Ra / R a c < 4 linear dynamics dominate two-dimensional flow behavior and the preferred wavelength decreases as Ra increases, but as R a / R a c increases above 4 the preferred wavelength increased (i.e. dA./d(Ra) < 0). For R / R c > 600 a preferred aspect ratio of 9.43 was obtained. Note that tests were conducted using an (eddy) Prandtl number of 1. (Note that the eddy Prandtl number ( G E ) represents the ratio of the eddy viscosity (v e ) to the eddy thermal diffusivity ( K e ) . The eddy Prandtl number will be discussed in more detail in sections 2.5.2 and 3.9 in association with the application of the model to the atmosphere.) However, numerical finite-difference simulations conducted by Sykes and Henn (1988) at high Rayleigh numbers did not result in an observed increase in scale. The discrepancy between these results and results of Rothermal and Agee (1986) was attributed by Sykes and Henn (1988) to sensitivity of the scale selection mechanism to the choice of numerical method, i.e. quasi-spectral versus finite-difference. Chapter 2: Background and Model Description 17 In summary, experiments indicate an increase in scale for values of the Rayleigh number above critical. However, both linear and weakly nonlinear theory are not able to account for these observations. Thus the question of the preferred scale at supercritical values of the Rayleigh number, in association with classical R-B convection, remains unresolved. Numerical modeling efforts have led to conflicting results as the ability of a two-dimensional nonlinear mechanism to account for observed increases in scales. We shall see in section 3.8.2 that for the case of classical R-B convection, current model simulations conducted at slightly supercritical values of the Rayleigh number for Prandtl number tests at 0.1, 1. and 10., indicate an increase in scales only for a = 0.1. These results are in agreement with those of Getling (1983) and Chang and Shirer (1984). Further investigations would be required in order to compare model results at high Rayleigh numbers with those of Rothermal and Agee (1986), and Sykes and Henn (1988). 2.3 Field Observations of Horizontal Roll Vortices The previous section outlined observations and theoretical developments which were initiated by controlled laboratory experiments of convection around the turn of the century (Benard, 1900), however evidence of geophysically occurring H R V accumulated more slowly. The pioneering observational work of Woodcock (1942) and Langmuir (1938) were instrumental in establishing the existence of coherent linear flow patterns in the form of longitudinal roll vortices within the atmospheric boundary layer and the oceanic boundary layer (where roll circulations aligned with the surface mean wind direction are denoted Langmuir circulations), respectively. The majority of field studies associated with H R V have been conducted over water which has the advantages of both surface homogeneity, and the verification of the existence of rolls from Chapter 2: Background and Model Description 18 satellite images of cloud street patterns, e.g. G A T E 1974 (LeMone and Meitin, 1984); KonTur 1981 (Briimmer, 1985); M A S E X 1983 (Atlas etal, 1986); M I Z E X 1983 and A C E 1984 (Hein and Brown, 1988); and ARKTIS 1988 (Briimmer et al, 1992). However, the occurrence of rolls over land is possibly more frequent as an increase in sensible heat flux which accompanies daytime solar heating may provide the necessary thermal instability, in the presence of sufficient mean wind shear, to initiate horizontal roll formation. However, due to complicated topographical effects and perhaps insufficient moisture to form clouds, clear-air roll formation over land has not received the same extent of field investigation as its over-the-water counterpart. Over-land studies include: the tetroon study of Angell et al (1968); M E T R O M E X (Kropfli and Kohn, 1978); JAWS (Christian and Wakimoto, 1989); as well as Reinking et al (1981). Over-land and over-water field studies have lead to the quantification of the following H R V characteristics presented in table 2.3 (from Etling and Brown, 1993). Property Observed Range vertical extent wavelength aspect ratio downstream extent orientation wrt geostrophic lifetime mean wind speed 10-1000 km -20 to +30 degrees 1-72 hours > 5 m/s 1-2 km 2-20 km 2-15 Table 2.3: Summary of observed roll characteristics. (Etling and Brown, 1993) Chapter 2: Background and Model Description 19 Other less conclusive aspects include the influence of topography, conditions favoring the formation of large aspect ratio rolls and multiple scales, the role of the clouds or cloud decks, and the degree of interaction with (and modification of) the mean wind. Since the current study focuses on scale selection, table 2.4 summarizes observed H R V aspect ratios which have been presented in the literature. Of particular interest are field observations where large aspect ratio H R V have been observed. These include H R V which are dominated by a single scale such as those of Briimmer (1985), Briimmer et al. (1992) and Kristovich (1993), as well as multiple-scale observations of LeMone and Meitin (1984), Walter and Overland (1984), and Hein and Brown (1988). Walter and Overland (1984) observed a number of scales of motion (1.3-1.7 km, 5-6 km, 12-15 km and 30 km) over the Bering Sea using aircraft, sonde, and satellite. The depth of the near-neutral boundary layer was approximately 0.7 km thus suggesting roll aspect ratios of ~2, 10, 20, and 40. Particular scales were found to dominate in different regions of the domain, and often smaller scales were superimposed on larger scales. The 0.3-0.5 km scales were attributed to turbulent plumes; 1.3-1.7 km scales to the inflection point instability; 5-6 km scales were classified as intermediate; 12-15 km scales were attributed to the influence of upstream topography; and finally 25-30 km scales were attributed to a possible resonant sub-harmonic of the smaller scales. Observations presented by LeMone and Meitin (1984) indicate the existence of cloud bands with horizontal wavelengths of 15-25 km. which were associated with a well-mixed shallow (-600 m) boundary layer suggesting aspect ratios of 30-50. However, unlike observations of Walter and Overland (1984) there was no apparent topographical source to account for the observed cloud spacings. Chapter 2: Background and Model Description 20 Reference A, (km) L / H Location Platform Comments Angell et al. (1968) 4 2.9 Idaho Ra buoyancy and shear Asai (1966) 20 7-10 Japan Sea Ai Atlas etal. (1983) 3-10 3-7 N.Y . and N.J. Sa C A O Atlas et al. (1986)1 1-2 1.3-3.3 East Coast of N.A. L i C A O Berger and Doviak 4 2.7 Norman, OK Ra (1979) Brurnmer (1985) 0.8-6.0 1.3-6.2 North Sea KonTur Brummer et al. 1.2-7 3.2-12.7 West of A i & Bo C A O (1992)2 Spitsbergen Christian and 5 avg 1.6 NE Colorado Ra JAWS; uni-Wakimoto (1989)3 directional wind Hein and Brown 2.3; 5.8 3.3; 8.3 Bering Sea Ai multiple (1988) scales; C A O Hildebrand (1980) 4.0-4.8 3.4-4 Central, O K Ra clear air rolls Holroyd(1971) 7-15 avg 5 Great Lakes Sa Kelly (1982) 3.5-7.5 2.7-5.8 Lake Michigan Ra snow storms Kelly (1984) 1.5-13.7 1-9.1 Lake Michigan Ra cold air flow Konrad (1968) 1.5 1.8 Wallops, V A Ra clear air;uni-directional wind Kristovich (1993) - 2.0-6.7 Lake Michigan R a & A i snow storms Kropfli and Kohn 7-9 4.5-6.5 St. Louis Ra urban (1978) environment Table 2.4: Summary of roll observations presented in the literature. Data collection methods: radar (RA), satellite (Sa), sonde (So), aircraft (Ai), tower (To), cloud observations (Clobs), smoke screen (Smoke), lidar (Li) and research vessel (Bo). Other abbreviations: cold air outbreak (CAO), Convection and Turbulence Experiment (Kontur), Joint Airport Weather Studies (JAWS), GARP Atlantic Tropical Experiment (GATE). 1 Roll wavelengths as observed by lidar. Nearby cloud streets had observed wavelengths that ranged from 5-10 km. 2 Downstream increase in aspect ratio (as observed by aircraft) from 3.2 to 5.7-12.7 over 200 km. 3 Predominantly unidirectional wind field within the PBL from 155°. Prediminantly unidirectional mean wind field above the inversion from 280°. Strong directional wind shear at the inversion with weak speed shear. Average PBL wind speed is 4.4 m/s which is much less than previous studies which range from 6 to 20m/s, and typically over 10 m/s (see table 1 of Christian and Wakimoto, 1989). Chapter 2: Background and Model Description 21 Reference A, (km) L / H Location Platform Comments Kuettner (1959) 7.5-10 2.5-4 Boston CLobs uni-directional wind Kuettner (1971) 2-8 2-4 ' Tropical Atlantic A i & S a uni-directional wind LeMone (1973) 1-7 2.1-6.5 Oklahoma & Great Lakes T o & S a slightly stable lapse rate LeMone and Pennell 2 2-3 Puerto Rico Ai marine layer (1976) LeMone and Meitin 15-25 25-50 East Atlantic A i multi-scales; (1984) Ocean G A T E Malkus and Reihl 4 2-3.5 Between Hawaii Ai (1964) & Philippines Markson (1975) 2-5 avg2 Maryland A i Meffiet al. (1985) 2-4 Atlantic Ocean Li C A O Miura(1986) 5.8-30 5.8-12.2 Eastern Asia Sa& So C A O Pitts et al. (1977) 1.5 1.5 Oklahoma Skylab Plank (1966) 0.3-5 avg2 Florida CLobs Puhakka and 2.5-10 3-10 Finland Ra Clear air Saarikivi (1986) rolls; rolls in rain & snow Reinking et al. 2-6 avg 2.7 Chickkasha, O K Ra (1981) Walter (1980) 3-12 2-8 Bering Sea Sa C A O Walter and 1.3-1.7; 1.8-2.4; Bering Sea A i , Sa C A O Overland (1984)4 12-15 17-21 Weston (1980) 2.5-9 2-5 England Sa Table 2.4 (cont.): Summary of roll observations presented in the literature. Data collection methods: radar (RA), satellite (Sa), sonde (So), aircraft (Ai), tower (To), cloud observations (Clobs), smoke screen (Smoke), lidar (Li) and research vessel (Bo). Other abbreviations: cold air outbreak (CAO), Convection and Turbulence Experiment (Kontur), Joint Airport Weather Studies (JAWS), GARP Atlantic Tropical Experiment (GATE). 4 Multiple scaling observed at wavelengths of 2 km (Ai), 12-15 km (Ai), and 30 km (Sa). Embedded in the 30 km bands of cloud streets were 5-6 km scales (Sa, Ai). Chapter 2: Background and Model Description 22 Observations of Hein and Brown (1988) also suggest the existence of multiple-scale H R V circulations. Roll aspect ratios were found to decrease with distance for the leading edge of an ice sheet while cloud street wavelengths increased. The cloud-street wavelengths were 2-4 times larger (satellite photos) than roll wavelengths calculated from spectra obtained using aircraft. Results of Hein and Brown suggest that the relationship between observed cloud street spacings and roll circulations may not be as direct as previously assumed. It is further suggested that cloud spacing resulting from an interaction of scales may account for these discrepancies. These observations are in contrast to observations of e.g. Miura (1986) and Brurnmer et al. (.1992) where increases in aspect ratio of H R V in the downstream direction were observed. Estimates of roll spacing using satellite photos has assumed that measurements of cloud spacing and an estimate of boundary layer depth using radiosondes is sufficient to determine the scale of HRV. However, observations of LeMone and Meitin (1984) and Hein and Brown (1988), suggest that cloud street scales may not (always) directly reflect the scale of the underlying H R V circulation. Christian and Wakimoto (1989) found that the clouds which formed in association with clear-air H R V , existed within the overlying stable layer, suggesting possible errors when estimating the depth of the boundary layer based on cloud top height. Also, the clouds were aligned with the mean P B L wind (unidirectional at 155°) rather than the mean stable layer wind field (unidirectional at 280°) within which they existed. Chapter 2: Background and Model Description 23 Gravity Wave Observations In addition to the large-scale and multiple-scale observations of H R V noted above, field observations relating to wave formation resulting from convective activity within the P B L are also of interest. Observations of Kuettner et al. (1987) (summarized in table 2.5) are presented in some detail as these observations and companion papers by Clark et al. (1986) and Hauf and Clark (1989) form the basis for the design of the model used in the current study. This over-land field study was designed to investigate internal gravity wave generation by the convective planetary boundary layer via a series of observational flights and numerical modeling efforts (Clark et ai, 1986; Hauf and Clark, 1989). It should be noted that cloud streets were observed only on 26/06/85 and 03/07/84, with vertical profiles of mean wind speed containing the typical low level jet within the PBL, reaching maximum velocities of 14 and 9 m/s, respectively. Vertical profiles of mean wind direction show a predominantly unidirectional mean wind field within the boundary layer, accompanied by a change of 70 and 125 degrees in mean wind direction across the inversion for 26/06/85 and 03/07/84, respectively. An aircraft investigation into the existing cumulus cloud on 03/07/84 found a change in wind direction from 280 to 50 degrees outside (above) to inside the cloud accompanied by a 5 m/s decrease in the wind speed. The resulting relative motion of the surrounding air against the cloud was from 265 degrees at 9 m/s. Thus Kuettner et al. suggested that the clouds may act as obstacles to the flow. Of particular note are the observations obtained on 12/06/84. Although this day was associated with an observed three-dimensional cumulus field, these atmospheric conditions form Chapter 2: Background and Model Description 24 the basis for the velocity and temperature profiles used in the two-dimensional numerical simulations of Clark et al. (1986). Date H i A0 i S z Ujnv Aw A c WM A X w d i r (km) (K) (gkg1) ( i o - V ) (m/s) (km) (km) (m/s) 12/06/84 2 3 4 7 10 5-6 10-11 3 E-SW 03/07/84 1.8 2.5 1.5 4 4 8 4; 9 2 50-275° 26/06/85 1.6 3 1.5 10 11 11-12 9-11 2 75° change Table 2.5: Observations of Kuettner et al. (1987). Abbreviations: H depth of the convectively driven PBL; A0 ; potential temperature discontinuity at the inversion; Aq; specific humidity discontinuity at the inversion; S zmean vertical shear; U i n v velocity at the inversion; A,w approximate horizontal wavelength of gravity waves; A,c approximate cloud spacing; w m a x maximum vertical velocity; W d i r wind direction within-above the boundary layer. However, the data set of Kuettner et al. (1987) is incomplete as aircraft observations of roll aspect ratio measured within the boundary layer are not presented. Also of interest would be measurements of gravity wavelengths measured higher up in the troposphere in order to compare with numerical results of Clark et al. (1986) and the current model. In conclusion, observations indicate the existence of rolls within a multitude of atmospheric environmental conditions ranging from: conditions of clear-air to snow storms; convectively unstable to near neutral environment; negligible cloud formation to complete coverage by a cloud deck; as well as roll formation over-land and over-water. It is therefore not surprising that a number of mechanisms have been suggested to account for H R V formation under specific conditions. However to date, no one explanation is able to account for all Chapter 2: Background and Model Description 25 observations of HRV. These H R V formation mechanisms, as well as roll modification mechanisms which have been proposed to explain large aspect ratio roll formation, are discussed in the next section. 2.4 Roll Formation and Modification Mechanisms How applicable are laboratory and theoretical results of section 2.2 to the atmosphere? This question may be addressed in part by considering the fact that H R V may exist in a convective atmosphere at Rayleigh numbers for which laboratory and theoretical studies of classical R-B convection indicate that only cellular planforms are stable. (Note that an observed range of Rayleigh numbers associated with the atmosphere are presented in section 3.9. Observations of Miura (1986) suggests that rolls may exist at Rayleigh numbers which are 2-3 orders of magnitude higher than laboratory observed values for stable roll formation.(Table 2.1)) This suggests that a stabilizing mechanism, in association with convection, must exist in order to produce stable roll circulations at these high values of the Rayleigh number. As early as 1942, Woodcock found that cellular convection was preferred when the mean wind speed was less than ~7 m/s and rolls when the mean wind speed was greater than ~7 m/s. These early observational results and others since (section 2.3) suggest that a mean wind speed in excess of ~5m/s is required in order to maintain stable roll circulations. Roll formation which is initiated and maintained by both buoyant and shear energy sources is termed shear-organized convection. Defining a right-handed coordinate system (x,y,z) with z positive upwards and y pointing in the downstream direction, plane parallel shear flow refers to flow associated with a mean velocity profile (u, v, w) of the form (0, v,0). Shear-organized convection (which incorporates the Chapter 2: Background and Model Description 26 effects of unstable thermal stratification on plane parallel shear flow) has been investigated for the case of constant mean wind shear ( d 2 v / d z 2 =0) (Kuo, 1963; Hill , 1968; Asai, 1970a) and variable mean wind shear ( d 2 v / d z 2 * 0 ) (Asai, 1970b; Kuettner, 1971). In particular Asai (1970b) used a linear perturbation analysis and found that for small values of the Richardson number Ri (where Ri = gccHAT / u *, and g is the acceleration due to gravity, a is the thermal expansion coefficient, H is the distance between the two plates, AT is the temperature difference between the upper and lower surface, and u* is a characteristic velocity scale for the base flow) an inertial instability associated with the form of the velocity profile dominates, while for larger values of Ri a three-dimensional longitudinal mode which is associated with the thermal instability dominates: i.e. rolls oriented at an angle to the mean wind direction in order that energy may be drawn from the mean flow. In addition to H R V which exist within a convectively unstable environment, rolls which are aligned at an angle to the mean wind may exist within a neutral environment (i.e. in the absence of buoyant effects) in association with an unstable mean velocity profile. As early as 1880, Rayleigh found that velocity profiles containing an inflection point were inherently unstable to infinitesimal perturbations. Investigations into atmospherically observed H R V have focused on the Ekman velocity profile (Ekman, 1905) which is the result of an (idealized) balance between the mean flow pressure gradient force, the Coriolis force, and viscous forces. The Ekman velocity profile includes speed and directional shear as both of the horizontal velocity components (u,v) are functions of height. The stability of the Ekman profile has been investigated both experimentally (e.g. Faller, 1965) and theoretically (e.g. Faller and Kaylor, 1966; Brown, 1970; Asai and Nakasuji, 1973). In general perturbation results indicate that the Ekman velocity profile Chapter 2: Background and Model Description 27 is unstable in the Reynolds number regime pertaining to geophysically observed flows (Brown, 1970). Incorporating unstable stratification into studies of the Ekman velocity profile, linear perturbation results of Brown (1972) indicate a shift towards smaller angles of orientation of the roll circulations with respect to the geostrophic wind direction, and a decrease in the preferred aspect ratio, with increasing thermal instability. (Note that the geostrophic wind results from the theoretical balance of the pressure gradient and Coriolis forces.) In general, both shear-organized convection and the inflection point instability associated with the Ekman velocity profile, predict aspect ratios which are of the order of 3 (Kelly, 1984: LeMone and Meitin, 1984) and are unable to account for observed large scale H R V . Thus, in order to explain large aspect ratio roll formation a number of roll modification mechanisms have been proposed including: anisotropic eddy viscosity (Priestly, 1962); large scale upwelling and subsidence (e.g. Krishnamurti, 1975 a,b,c); boundary conditions (e.g. Sasaki, 1970); latent heat release (e.g. Sykes et al., 1988); wave-wave interaction (Mourad and Brown, 1990); and gravity wave-boundary layer interaction (e.g. Clark et al., 1986). These mechanism are discussed below. Anisotropy of Eddy Viscosities Priestly (1962) attributed discrepancies between observed cellular aspect ratios (-30) and linear theory (~3) to anisotropy of sub-cell scales of turbulence. In particular it was assumed that the horizontal (eddy) viscosities are an order of magnitude or two greater than the vertical (eddy) viscosity. Priestly notes that such inequalities are not true of general turbulence but says that eddy viscosity anisotropy may be a property under stable conditions. Chapter 2: Background and Model Description 28 Miura (1986) incorporated the effects of anisotropic eddy viscosity proposed by Priestly (1962) into the linear perturbation model of Asai (1968) for cellular roll convection. It was found that viscosity ratios of 10-20 were required in order to produce roll aspect ratios of 5-10. However, there is no observational evidence to support the concept of anisotropic eddy viscosities in boundary layers with neutral or unstable thermal stratification (LeMone and Meitin, 1984). Large Scale Subsidence / Upwelling A three part series by Krishnamurti (1975 a,b,c) investigated the relationship of open or closed cellular cloud patterns with respect to large scale sinking (W < 0) or rising (W > 0) motion (where W is the vertical component of the mean wind). A weakly nonlinear model was developed (1975a) and results indicate that rolls would be stable for a sufficiently large lapse rate and sufficiently small vertical velocity, while hexagons would be stable for W sufficiently large and a sufficiently small lapse rate. A laboratory model was developed (1975b), and the regimes for which these two modes of convection were theoretically predicted to be stable, were found to be in good agreement with experimental results. In agreement with perturbation and lab results of Krishnamurti (1975 a,b,c), three-dimensional finite-difference model results of Chlond (1992) also suggest that observed values of large-scale subsidence are too small to have a significant effect on roll dynamics. Chapter 2: Background and Model Description 29 Thermal Conductivity of the Boundaries and Boundary Conditions Application of constant temperature boundary conditions, assumes that the thermal conductivity of the bounding surfaces ( K S ) is infinite or at least much greater than the thermal conductivity of the fluid ( K F ) . Linear perturbation investigations into the effect of finite conductivity of the bounding planes on R-B convection indicate a significant decrease in R a c and an increase in X, c , as K F / K S increases (Hurle et al., 1966; Nield, 1968; Proctor, 1981). In particular, Hurle et al. (1966) found that for K F / K s = [0,°°) the value of the critical Rayleigh number for the case of no-slip boundary conditions was found to decrease from -1708 to 720. The corresponding preferred wavelengths were 2.016 increasing 3.033 for K F / K S = 2 and to infinity as K F / K S - > < » . Sasaki (1970) used boundary conditions which he felt were more representative of conditions relating to the atmospheric boundary layer and incorporated thermal influences associated with the surface layer. Sasaki investigated large aspect ratio of cellular convection using surface and inversion boundary conditions which were expressed as a linear combination of the vertical temperature gradient and the temperature itself. Results were found to be largely dependent on the ratio of the vertical temperature gradient to the temperature of the boundaries. As with the studies relating to R-B convection, a decrease in R a c and an increase in A,c was observed as the ratio of the radiative heat flux to the conductive heat flux at the boundaries decreases. Although the surface boundary conditions may be more 'realistic', application of these same conditions to the inversion is questionable. Chapter 2: Background and Model Description 30 Latent Heat The effects of latent heat release on boundary layer dynamics has lead to conflicting results with respect to large aspect ratio roll formation: Two-dimensional numerical results from the finite-difference model used by Mason (1985), suggest that aside from the importance of cloud-top radiative cooling which prohibits rapid cloud deck formation, boundary layer dynamics incorporating latent heat effects were similar to that for a dry C B L . However, two-dimensional finite-difference model results of Sykes et al. (1988) indicate that aspect ratios of H R V larger than 4 were not possible without incorporation of latent heat effects and in particular, cloud top entrainment. Under certain conditions aspect ratios of up to 10 were successfully modeled. However, they found that large scale organization did not occur until clouds were deep enough in order that latent heat release became significant and entrainment fluxes dominated over surfaces fluxes. Their model was assumed applicable for H R V formation over water which is in general associated with a larger latent heat flux and (possibly) smaller surface heat flux, when compared to roll formation over land. Chlond (1988) obtained roll aspect ratios of 3-18 using a two-dimensional finite-difference model. The model was based on ascending motions which are dry adiabatic below the cloud condensation level, moist adiabatic above, and all descending motions were dry adiabatic. For an eddy Prandtl number of 0.5, it was found that the aspect ratio of boundary layer rolls increased with increasing static stability. This result agrees well with observations which indicate an increase in H R V aspect ratio and an increase in stability in the downstream direction (e.g. Miura, 1986). Chapter 2: Background and Model Description 31 An extension of the two-dimensional, seven-coefficient spectral model of Chang and Shirer (1983) (section 2.2) was used by Laufersweiler and Shirer (1989) to investigate latent heat effects for the case of a stratocumulus topped boundary layer. Al l motion above the base of the cloud deck was assumed moist adiabatic, and all motion below the base of the cloud deck was assumed dry adiabatic. For the case of ~ 10% cloud cover, results indicate that the effects of latent heat release is small. Roll circulations were symmetric and the preferred aspect ratio was similar to that obtained for dry convection. For larger cloud decks (1/2 to 1/3 of the domain) the flow circulation within the cloud layer was found to detach from the sub-cloud layer. The circulations were found to become increasingly asymmetric as latent heating increased. Results of these studies into the effect of latent heat release on roll dynamics suggest that such effects are not able to account for large aspect ratio roll development under conditions of negligible cloud formation. Wave-wave Interactions An explanation of multiple-scale roll development resulting from triadic wave-wave interaction in a neutral environment was presented by Mourad and Brown (1990) using a nonlinear spectral model. Buoyant effects were neglected and roll formation was the result of the inflection point instability associated with a mean Ekman velocity profile. Both wave-wave and wave-mean flow interactions were modeled in an attempt to explain large aspect ratio roll development. Model results were compared to observations of Walter and Overland (1984). The analysis by Mourad and Brown includes three principal waves and the first harmonics which were assumed to be forced by the principal wave and were not free to develop on their own. The harmonics were Chapter 2: Background and Model Description 32 considered to be important for the energetics of the interactions but contributed little to the spatial structure of the large eddy field. Mourad and Brown suggest that the determining factor in the selection of a large eddy state is the presence or absence of source waves (which may result from other instability mechanisms or topography) that would project onto the eigenstates of the boundary layer, producing well defined scales in the flow field. However, the restriction on the interaction between the harmonics and other scales limits natural flow development. Results from the current model presented in chapter 3, suggest that in fact the harmonics of the principal waves play a significant role when full interaction between scales is allowed. Gravity Wave-Boundary Layer Interactions Gravity waves were generated in the stable upper region above the convective boundary layer in the two-dimensional finite-difference model of Mason and Sykes (1982) which was designed to investigate the unstable PBL. However, no influence of the gravity waves on boundary layer dynamics was observed. In contrast, two-dimensional numerical results from the finite difference model used by Clark et al. (1986) indicate a strong coupling of the stable upper layer and boundary layer dynamics. In particular it was found that the result of the coupling was a feedback of the preferred (longer wavelength) gravity wave scale on the scale of convection within the boundary layer; i.e. the gravity waves 'tuned' their own energy source. The result was an increase in observed H R V scales as compared to those for a boundary layer which is considered in isolation. An extension to three-dimensions by Hauf and Clark (1989) in general corroborated the two-dimensional results. Chapter 2: Background and Model Description 33 Results of Hauf and Clark also indicate that the difference in preferred angles of orientation between the rolls within the boundary layer and the overlying gravity waves results in an overall broken pattern. The apparent absence of significant interaction between the induced gravity wave field with boundary layer dynamics as observed by Mason and Sykes (1982) has been suggested by Clark et al. (1986) to be the result of weak shear (~ 1.4xl0~3 / s) and weak surface heat flux values (~ 30watts / m 2 ) used in their model simulations. Sang (1991, 1993) developed an analytical two-layer model using linear perturbation techniques and found that gravity waves induced by a mechanical or thermal disturbance were composed of trapped waves with selected wavenumbers. These wavenumbers were determined by the interaction between the convective boundary layer and the overlying stable layer. The explicit relation of wavenumber selection defines the pattern of the vortices which may appear as cloud streets, cellular convection or convection thermals. However, results from the current study will show (section 3.8) that a linear analysis biases results and may lead to incorrect scale selection. Balaji et al. (1993) used numerical finite-difference simulations to investigate tropical cloud clusters as observed during the 1974 G A T E field program (LeMone and Meitin, 1984). Balaji et al. (1993) attributed observed cloud band widths (of ~ 25 km) to the interaction of three distinct modes: the Rayleigh mode of forced convection within the PBL; a nonhydrostatic trapped mode in the lower atmosphere; and a freely propagating mesoscale wave mode existing within the entire troposphere. It was further suggested that the subsidence associated with the large scale (deep) gravity waves prohibited cloud formation within regions of (large scale) descent. However within the ascent regions, enhanced cloud formation results from constructive interference of the three dominant scales. Chapter 2: Background and Model Description 34 2.5 Physical Model Description As outlined in section 2.4, a number of mechanisms have been presented in the literature which attempt to explain large aspect ratio modes of convection however, no one mechanism is able to account for all conditions under which H R V have been observed (section 2.3). Motivated by observations of Kuettner et al. (1987) and numerical modeling results of Clark et al. (1986), the current model investigates coupled interactions of large scale gravity waves and convective boundary layer dynamics as a possible mechanism for the formation of large aspect ratio rolls. Apart from large scale-subsidence (which has been shown to be unimportant (section 2.4)) all other mechanism discussed in section 2.4 may be considered as local or internal boundary layer process. In contrast, the hypothesis that the scale of convection may be influence by induced tropospheric gravity waves suggests the importance of large scale external boundary layer processes and the need to consider the full depth of the troposphere when modeling the PBL. In contrast to the numerical modeling efforts of Clark et al. (1986) and Hauf and Clark (1989), the current study adopts a more theoretical approach which attempts to isolate possible gravity wave influences on the convectively driven PBL by removing any complicating effects such as dynamic instabilities, latent heat release, thermal conductivity properties of the boundaries and entrainment. Referring to figure 2.1, it has been proposed (Clark et al., 1986) that the updrafts/downdrafts associated with H R V warp the inversion base. The warped interface then acts as a series of obstacles to the flow above the boundary layer and internal gravity waves are generated in the overlying stable layer at a wavelength determined by the convective activity. It is Chapter 2: Background and Model Description 35 further suggested that the induced gravity waves are influenced by tropospheric characteristics (such as stability and mean wind) and as a result, the wavelength of the gravity waves increases. This larger horizontal scale is then fed back through mechanical/thermal forcing at the inversion, thus influencing the scale of convection within the boundary layer. Based on this hypothesis, the current physical model consists of a two-layer model of the lower atmosphere (following that proposed by Clark et al., 1986) which extends from the ground to a height just above the tropopause (approximately 10-15 km). The lower layer of the physical model represents the fully developed PBL which forms the lowest 1-2 kilometers of the Earth's atmosphere, and is driven by daytime surface heating. The upper layer of the physical model represents the remainder of the troposphere which is characterized by increasing potential temperature with height and is therefore stable. The upper layer extends from the top of the PBL up to the tropopause and into the lower portion of the stratosphere. Although the importance of the dynamic coupling of the convective boundary layer and the overlying stable layer has been noted, the current theoretical modeling approach involves considering each of these two regions independently but incorporating dynamic and/or thermal influences of one regime into the other. Before presenting a detailed description of the two model domains, the assumptions used to develop the current model and its regime of validity are presented. Chapter 2: Background and Model Description 36 TROPOPAUSE INVERSION Figure 2.1: Physical model of the troposphere. Highly idealized as the amplitude of the gravity waves may vary with height. Note that the wave field has not yet propagated to the tropopause. (See text for model description.) Chapter 2: Background and Model Description 37 2.5.1 Assumptions and Regime of Validity In order to develop the current model, the following assumptions are applied: • Curvature effects of the Earth are neglected. Additional terms resulting from the expansion of the equations in spherical coordinates are omitted and the local (x,y,z) coordinates replace the spherical coordinates (cp,0,r). • The only external force applied is that due to gravity. Therefore F is replaced by (-gk) in (2.2). • The Boussinesq approximation is applied (Boussinesq, 1903). Thus i) the variation of perturbation density is neglected in the conservation of mass equation as well as in the momentum equation except in association with buoyancy effects; and ii) the influence of pressure perturbations on buoyancy are neglected. See Mahrt (1986) for conditions under which application of the Boussinesq approximations is considered. ° Thus the equation of state (2.4) is written in the form p = p 0 ( l - c c ( T - T 0 ) ) (2.9) where a is the coefficient of thermal expansion and the subscript 'o' denotes evaluation at the reference temperature; p is replaced with p 0 in association with the pressure terms in the momentum equation and by (2.5) when p is multiplied by gravity. 0 Also the fluid is assumed incompressible (i.e. dp / dt = 0) and (2.1) reduces to V V = 0 (2.10) • The total pressure (P) is assumed to be a linear combination of a constant (p 0) and a variable pressure (p) of the form P = p 0 + p. The constant pressure represents a hydrostatic Chapter 2: Background and Model Description 38 balance of the atmosphere which may be obtained using (2.2) and (2.10), and is given by dp 0 /dz = - p 0 g . • Based on the results of Laufersweiler and Shirer (1989) and the desire to isolate gravity wave influences alone, latent heat effects are neglected. • Inertial effects due to the Earth's rotation are neglected. Thus eliminating flow instabilities associated with a mean velocity profile which turns with height. • The boundary layer is assumed to be unstable and the mean wind within the boundary layer is assumed to be unidirectional (e.g. Kuettner.et al., 1987). Two-dimensionality is imposed by further assuming that variations in the along axis direction are prohibited (i.e. 3 / dy — 0 where the coordinates are aligned such that y is oriented along the roll axis), resulting in the decoupling of the governing equations in the along wind and cross wind directions. Therefore the net assumption is that the flow (secondary circulation) derives all of its energy from buoyancy effects, and that there is no energy supplied by the mean wind. • In general, the C B L is characterized by a superadiabatic surface layer (bottom 5-10%), a well-mixed layer (middle 35-80%), and an entrainment zone (top 10-60%). (See Stull (1988) for characteristics associated with each of these three layers). The current study makes no attempt to model explicitly either the surface layer or the entrainment zone, and many interesting and challenging questions associated with both of these regions are actively under investigation. Thus only bulk characteristics of the mixed layer are considered. Since the non-local transport by large eddies is explicitly modeled, unresolved scales of turbulence are assumed adequately represented by a first-order closure method (K-theory), i.e. it is assumed that the effects of turbulence on the flow may be represented by an 'eddy viscosity' approach. White (1988) investigated Chapter 2: Background and Model Description 39 experimentally the effects of viscosity on flow development in association with R-B convection. Results indicate that the effect of a temperature dependent (or similarly a spatially dependent) viscosity was to reduce the effective depth of the fluid as a cool, viscous layer developed close to the upper boundary. The flow was observed to develop in 'less' viscous regions which resulted in a reduction of the preferred scale. These results suggest that the effect of a more 'realistic' eddy viscosity profile which has been proposed to attain a maximum value in the lower portion of the boundary layer (see O'Brien, 1970), may be to separate the large eddy circulation from the lower surface as flow development would be inhibited in this region. • Motivated by the resulting analytical simplicity, free-slip boundary conditions are applied. Although the free-slip boundary conditions do not accurately represent (at least) the bottom surface of the atmosphere (i.e. the ground, which would be better modeled using no-slip conditions), it was felt that the mathematical simplifications would still lead to valid results which may be used to further understand the dynamics associated with the physical problem. Results from linear theory of R-B convection (section 2.2) indicate that the effect of imposing no-slip boundary conditions is to reduce the effective depth of the layer which results in a decrease in the wavelength of the preferred scale at the critical Rayleigh number. Thus it may be anticipated that with respect to the current model, application of no-slip conditions at the surface would lead to a flow circulation which is detatched from the boundaries and exists within a fluid layer whose depth is (marginally) less than that of the boundary layer. • The upper and lower boundaries are assumed rigid. Undoubtedly satisfied at the lower boundary (i.e. the Earth's surface), assuming that the inversion may be represented by a 'rigid' boundary is certainly an idealization of the physical model as the inversion zone is unarguably Chapter 2: Background and Model Description 40 deformable. Since the effect of the growth and decay of the convective boundary layer is beyond the scope of the current study, it is assumed that the height of the mixed layer is relatively fixed. Then it may not be unrealistic to assume that fluid within the convectively driven boundary layer remains confined between the ground and the depth of the mixed layer, and that both the wind shear and atmospheric stability of the overlying stable layer prohibit significant mixing. Under these assumptions a sinusoidal forcing at the inversion may be interpreted as representing a traveling wave which is moving with a constant velocity equal to the constant mean wind field, at a time scale less than that associated with entrainment (or detrainment), and that the frame of reference corresponds to the reference frame of the traveling wave. • In order to eliminate the effects of imperfectly conducting bounding surfaces both the surface of the earth and the inversion are assumed to be at a constant temperature. Assuming that the inversion may be modeled as a constant temperature surface (or equivalently as a surface of constant density) is perhaps the most controversial of all of the applied boundary conditions, and is unarguably highly idealized. This boundary condition plays a leading role in the form of the solutions obtained, since the assumption that the upper surface is of a uniform temperature results in a periodic warping of the mean temperature field. Although it would be possible to devise a laboratory experiment which would simulate these conditions, it is questionable as to whether or not this is truly reproduced in the atmosphere. However, under the 'traveling wave' interpretation, the inversion is assumed to exist as a density discontinuity interface within the fluid, along which the gravity wave propagates. Assuming horizontal uniformity of the PBL depth (in the unperturbed state), the deformation of the inversion may then be envisioned as a surface of constant density and thus of constant temperature. Chapter 2: Background and Model Description 41 • In order to avoid the inadequacies of linear theory, a nonlinear approach has been adopted. Furthermore, to ensure full dynamical freedom of scale selection, a Fourier integral analysis of the horizontal flow dependence will be applied. • The mean wind above the boundary layer is also assumed to be unidirectional and consists of the component of the actual three-dimensional velocity field which is normal to the roll axis. Observational support for a (near) 90 degree shift in the mean wind direction above the boundary layer includes Jaeckisch (1972). • With respect to modeling of the stable upper layer, variations in the mean density with height are neglected. This assumption may be justified by noting that with respect to the dynamics of the upper layer, the current study is primarily interested in determining the dominant scales in the region of the inversion. Results of Wurtele et al. (1987), which include a background density profile which decreases exponentially with height, noted an additional term in the governing equations which filtered-out waves with wavelengths greater than -100 km. Thus the proposed model is similar to shear-organized convection models at high Richardson numbers, but with the added restriction of no along-axis variation of the flow variables, i.e. it is assumed that the roll axis is parallel to the mean wind direction. The absence of latent heat effects and an unstable boundary layer suggests applicability to H R V formed by daytime surface heating over horizontally homogeneous terrain with little or no cloud formation. Since the current model neglects the growth of the boundary layer and entrainment at the inversion, this implies that the current model would be most valid in association with the early afternoon hours when the depth of the boundary layer is least variable, i.e. before the afternoon Chapter 2: Background and Model Description 42 decay which accompanies the decrease in solar radiation and therefore a decrease in surface heat flux. It is also assumed that there do not exist alternate sources of gravity wave formation such as nearby topography (although the results of the current model apply to gravity wave coupling at the inversion regardless of the source of the gravity waves. The imposed two-dimensionality does however suggest that the axis of the wave field is parallel to the roll axis direction). 2.5.2 The Convectively Driven Planetary Boundary Layer In order to simulate the effects of imposing a larger scale at the inversion base, the lower layer is modeled by analogy with a modified R-B convection model in which the upper bounding surface is replaced with a corrugated lid of prescribed wavelength and amplitude as illustrated in figure 2.2. The imposed wavenumber is assumed to correspond to the preferred scale associated with the overlying gravity wave field, and both the wavelength and amplitude of the warping are assumed known. > C B L Ground Figure 2.2: The convectively driven planetary boundary layer. The effects of the gravity waves on boundary layer dynamics are modeled as a sinusoidal warping of the inversion at a wavelength A with an amplitude of e. Chapter 2: Background and Model Description 43 Analysis Methodology As noted in section 2.5.1 a nonlinear approach has been adopted. The inclusion of the advection terms CV • V)V and (V-V)T (where V represents the velocity vector with components (u,v,w), and T is the temperature) into the momentum equations and energy equation respectively, allows the transfer of energy between various scales of motion. Thus the nonlinear terms allow the flow to develop 'freely' at wavelengths other than that initially predicted by linear theory via the growth and decay of modes as energy is either transferred to, or removed from a particular scale. In contrast, linear wave theory only allows constructive and/or destructive interference resulting from a superposition of the wave field. The choice of analytical development versus numerically solving the original equations requires careful consideration of the objects of the study. The nonlinearity of the modified R-B convection problem will naturally lead to the use of numerical algorithms even if analytical techniques are initially applied. Thus the advantages and disadvantages must be carefully considered. Spectral methods and quasi-spectral methods (section 2.1) require initial analytical development of the equations and have been used to study classical R-B convection by (for example) Getling (1983), Chang and Shirer (1984) and Rothermal and Agee (1986). Numerical investigations into solutions of (2.1)-(2.4) have been conducted by a number of researchers using finite-difference schemes (e.g. Sykes and Henn, 1988) and results from these studies have been presented in section 2.2. Since the current study focuses specifically on scale selection, representation of the governing equations in a spectral or pseudo-spectral mode is naturally suggested. Spectral methods (applied to R-B convection) typically assume that the horizontal flow dependence may Chapter 2: Background and Model Description 44 be represented as a Fourier series (see section 2.1 and equation (2.7)). However, the concept of a fundamental mode and its harmonics arises from the dynamics of a linear system and the assumption that the nonlinear system would chose a length scale corresponding to an integral multiple of the preferred wavelength selected at the critical Rayleigh number has been suggested to be too restrictive (Getling, 1983). Thus the validity of imposing such a prescription on the allowable scales of nonlinear convection is questionable. The analytical techniques used in the analysis of the nonlinear flow behavior associated with the modified R-B convection model are based on the Fourier integral techniques used by Getling (1983). The use of the more general Fourier analysis is considered preferable to a series representation since nonlinear interaction of scales over a continuum in wavenumber space does not restrict wavenumber selection to that of the fundamental wavelength or its harmonics. Although vertical resolution is limited by the number of vertical modes included in the Fourier expansion of the flow variables, thus perhaps suggesting a pseudo-spectral analysis, careful consideration to ensure that the solution obtained is (relatively) insensitive to an increase in the vertical modes, allows confidence that the solution obtained captures the essential flow features. (Note that for high Rayleigh number flow (i.e. Ra > 400Ra c) the quasi-spectral model of Rothermal and Agee (1986) required 48 grid intervals in the vertical, and the finite-difference model of Sykes and Henn (1988) used 51, suggesting that a large number of vertical modes is necessary to resolve (in particular) the thermal structure.) Chapter 2: Background and Model Description 45 2.5.3 Stable Upper Layer Referring to figure 2.3, the effects of boundary layer dynamics on the overlying stable layer are incorporated by assuming that the effect of the periodic initial warping of the inversion by the boundary layer roll vortices is analogous to flow over topography. Since it is the perpendicular component of the mean flow which results in the generation of a wave field, this effect will be maximized when the wind field above the C B L is at right angles to the roll axes, i.e. perpendicular to predominant mean wind direction within the boundary layer. The stable domain is modeled using both a one-layer (section 4.4) and a three-layer (section 4.5) model, the latter incorporating more features of the assumed stability and velocity profiles. mean wind X Figure 2.3: The stable upper layer. Initial warping of the inversion (z = h) by the convectively driven roll vortices will occur at the wavenumber X. (Note that the vertical arrows denote updrafts and downdrafts associated with the roll circulation, and that the mean wind within the boundary layer is assumed into the page.) Chapter 2: Background and Model Description 46 Analysis Methodology A linear perturbation analysis of (2.1)-(2.4) as discussed in section 2.1, is applied to governing equations associated with model of the stable upper layer. The investigation of the response of the upper layer to topographic effects focuses on a particular atmospheric environment scenario. Specifically, the dynamics of the stable upper layer are investigated using stability and mean velocity profiles estimated from temperature and velocity profiles of Clark et al. (1986) which they used in two-dimensional numerical simulations (section 4.2). Of particular interest is the resonant response of the atmosphere and it's preferred wavelength of oscillation as well as the primary scale of influence at the top of the convectively active boundary layer. Details of the analytical development are presented in chapter 4. 2.5.4 Summary of Research Objectives The overall objective of the current study is to investigate the coupled interactions of the convectively driven PBL and an induced gravity wave field in the overlying stable layer as a possible mechanism for the formation of large aspect ratio (i.e. > 4) H R V . This will be done by meeting a series of sub-objectives: • To investigate the nonlinear flow dynamics associated with a modified Rayleigh-Benard convection model in which the flat upper bounding plane is replaced with a corrugated surface. The wavelength and amplitude of the warping are assumed known.; ° To develop analytical forms for the two-dimensional nonlinear governing equations using a weakly nonlinear perturbation analysis as well as a fully nonlinear analysis; Chapter 2: Background and Model Description 47 ° To transform the governing equations using Fourier integral techniques to wavenumber space; 0 To investigate the linear response of the systems of equations; ° To obtain numerical solutions to the nonlinear equations governing the Fourier coefficients; ° To apply the results of the modified R-B convection model to the atmosphere using atmospherically observed values of the parameters. • To investigate the linear response of the stable upper layer using a one-layer model. ° To determine the resonant mode of the system based on atmospheric profiles used by Clark et al. (1986), than is used in the one-layer model. ° To determine the dominant mode using kinetic energy scaling considerations. • To investigate the linear response of the stable upper layer using a three-layer model. ° To incorporate more vertical features of the velocity and buoyancy frequency atmospheric profiles of Clark et al. (1986); ° To determine the resonant modes of the three-layer system; ° To compare with results obtained from the one-layer model. • To examine techniques which may be used to dynamically couple the convectively driven planetary boundary layer and overlying stable layer. Chapter 3 Modified Rayleigh-Benard Convection and the Convectively Driven Planetary Boundary Layer 3.1 Modified Rayleigh-Benard Convection Model As discussed in section 2.5.2, the effects of a imposing a length scale at the inversion on convection within the PBL is modeled by analogy to a modified R-B convection model. The classical problem of thermally driven flow between two flat plates of infinite horizontal extent is replaced by a model consisting of a flat lower surface and corrugated upper surface, as depicted in figure 3.1. Modification of the classical problem involves replacing the upper bounding plane z = H by a corrugated surface defined by z = H - 1 + ecos (II , where H is the mean distance separating the two plates, 8 is the amplitude of the forcing at the upper plate (assumed small), and b is the imposed nondimensional wavenumber as selected by the induced tropospheric gravity waves, with e and b assumed to be known. The temperature of the lower and upper surfaces are denoted T 0 , and T i respectively, and in order that the convective instability may be initiated, it is assumed that the vertical temperature gradient is negative. Other a-typical boundary studies Although surface nonuniformities are an accepted practical limitation of experimental design, initial interest into possible advantages of a-typical geometries arose in association with 48 Chapter 3: Modified Rayleigh-Benard Convection ... 49 engineering flows; in particular their effect on heat transfer and the possibility of maximizing heat transfer efficiency. Watson and Poots (1971) studied the effects of sinusoidal, isothermal bounding surfaces on flow between vertical walls. In their model both endwalls are corrugated and 180 degrees out of phase. The current study follows the analytical development by Watson and Poots (1971); i.e. the choice of nondimensional variables; mapping of the domain to the 'classical' domain; and a weakly nonlinear perturbation expansion of the flow variables. 27tH/b Figure 3.1: Modified Rayleigh-Benard convection model domain. The classical domain is modified to include a warped, isothermal rigid upper lid. The mean distance between the two bounding surfaces is denoted by H . The amplitude of the warping (e) and the imposed wavenumber (b) are assumed known. Chapter 3: Modified Rayleigh-Benard Convection ... 50 Kelly and Pal (1978) investigate the effects of spatially periodic boundary conditions of flow between horizontal bounding planes using a weakly nonlinear perturbation analysis combined with a truncated spectral representation of the flow variables. They consider two sets of boundary conditions: flat plates with sinusoidal temperature; and warped plates with constant temperature. The amplitude of the sinusoidal variations (either thermal or spatial) are assumed small, however the wavenumber associated with the imposed boundary conditions is assumed to be identical to the critical wavelength associated with the initiation of convection of the classical R-B convection problem. A possible phase difference between the boundary conditions at the surfaces is included in the analysis. They found that a circulation of order e is induced for all values of the Rayleigh number less than critical. However for Rayleigh numbers near critical, consideration of the nonlinear terms resulted in an amplitude of the secondary circulation of order £ 1 / 3 and not unbounded as predicted by linear theory in association with forcing at the resonant wavelength. Although the current model is not unique in its geometry, it does not impose restrictions on the choice of the wavelength associated with the surface corrugation (c.f. Kelly and Pal, 1978). Analysis of the governing equations using both a weakly nonlinear perturbation analysis and a fully nonlinear approach will incorporate a more general form for the expansion of the flow variables and therefore considers a wider range of nonlinear interactions using Fourier integral techniques, when compared to the highly truncated expansion for the flow variables used by both Watson and Poots (1971) and Kelly and Pal (1978). The remainder of the chapter outlines the analytical and numerical techniques used to investigate the flow behavior associated with the modified R-B convection model. Chapter 3: Modified Rayleigh-Benard Convection ... 51 After applying some simplifying assumptions (section 3.2), the two-dimensionality of the problem is exploited via the introduction of a streamfunction, and the number of equations are reduced from six to two. Consideration of important length, time and velocity scales leads to the introduction of the Prandtl and Grashof numbers as the equations are reduced to their nondimensional form (section 3.2.2). In order to investigate the behavior of the solution for the streamfunction and temperature, the governing equations (section 3.2.2) undergo a series of transformations. The first involves a coordinate transformation (section 3.3) designed to simplify application of the boundary conditions (of section 3.4), followed by a Fourier integral and finite Fourier series analysis (section 3.5) which transforms the equations to wavenumber space. A linearized version of the resulting integrodifferential equations for the Fourier coefficients is investigated analytically (section 3.6). A discussion of the numerical solution procedure used to solve the complete, nonlinear integrodifferential equations (section 3.7) is followed by a presentation of the numerical results (section 3.8). Finally, in order to compare with atmospherically observed phenomenon, a numerical test was conducted using atmospherically derived values of the parameters, and is presented in section 3.9. The analysis procedures are summarized in figure 3.2 with the corresponding section numbers given in brackets. Chapter 3: Modified Rayleigh-Benard Convection . Governing Equations Physical Space (3.2) Nondimensionalized (3.2.1) Transformed to Mapped Space (3.3) Weakly Nonlinear Perturbation Analysis (3.3.2) Zeroth Order Equations (3.3.3) Fourier Transformed to Wavenumber Space (3.5.1) Numerical Solution (3.8.2) First Order Equations (3.3.4) Fourier Transformed to Wavenumber Space (3.5.2) Linear Analysis (3.6.2) Numerical Solution (3.8.3) I I Results of the Zeroth Order Equations and First Order Equations are combined to determine the net response as predicted by the Perturbation Analysis (3.8.5) Fully Nonlinear Truncated Analysis (3.3.2) Truncated Equations (3.3.5) Fourier Transformed to Wavenumber Space (3.5.3) Numerical Solution (3.8.4) I Response predicted by the Truncated Analysis (3.8.5) Figure 3.2: Flow chart outlining the analysis procedure. Chapter 3: Modified Rayleigh-Benard Convection ... 53 3.2 Governing Equations in Physical Space The governing equations presented in section 2.1 are further developed and adapted for study of the modified R-B convection model. In addition to the assumptions of section 2.5.2, it is also assumed that • there exists no mean flow in the cross-roll plane (by definition of R-B convection); • the temperature difference (T - T m ) is small compared to the mean temperature (Tm) between the upper and lower surfaces where T m = (T0 + T,) / 2; • physical parameters such as the thermal diffusivity (K), the kinematic viscosity (v), etc. are assumed constant over the temperature range considered; • and the equation of state (2.4) (or 2.9) may be written in the form ( P ~ P m ) = ( T - T J p T 1 ; where the subscript'm' denotes evaluated at the mean temperature. Under these assumptions the (2.1) through (2.3) reduce to — + ( v v ) v = - — V p + v V 2 V + g ( T ~ T m ) k (3.2) dt V ' p m T m r m m — + ( V - V ) T = K V 2 T (3.3) at v ; V - V = 0 (3.4) where we have used (3.1) to eliminate density from (3.2). Note that V = (u,w) where u and w represent the horizontal and vertical component of velocity respectively. Also, K = K M / (pc ) is the thermal diffusivity, and both K and v are evaluated at the reference temperature (Tm), p is the Chapter 3: Modified Rayleigh-Benard Convection ... 54 dynamic pressure (i.e. the difference between the actual pressure and the hydrostatic pressure in the absence of heating), and k is the unit vector in the z-direction. Since we are considering a two-dimensional analysis, the equation of continuity (3.3) suggests the introduction of the scalar streamfunction which is defined in relation to the horizontal and vertical components of velocity as u = - ^ a n d w = - ^ . (3.5) dz dx Eliminating the pressure terms from the momentum equations, the equations governing the velocity components and temperature reduce to the following for the streamfunction and temperature V T m J T.V72 \ V74 3 T T - T . ^ — J(V>,\ | / ) = v V > - g — dt dx 3T (3.6) ^ - J ( X | / , T ) = K V 2 T . (3.7) dt where J = [f pg q - f q g p j is the Jacobian and the subscripts denote differentiation with respect to the coordinates (p,q) where for (3.6) and (3.7) (p,q) represent (x,z). Equations (3.6) and (3.7) may be written in an alternate form by introducing vorticity defined by Z, = V x V , which is expressed in component form as ( £ , , £ 2 , £ 3 ) = [ w y - v z , u z - w x , v x - u y ] , where again the subscripts denote partial differentiation with respect to x, y, or z. Therefore the term V 2 \ | / is equivalent to C,2 and (3.6) and (3.7) may be written Chapter 3: Modified Rayleigh-Benard Convection 55 a t + V - V ^ 2 = v V 2 C 2 - g T ^ dx a ( T - T . A V Tm J — + V ' V T = K V 2 T . a t (3.8) (3.9) In this form it may be readily seen that the terms from left to right in (3.6) (or 3.8), represent the local rate of change of vorticity, advection of vorticity by the velocity field, viscous dissipation of vorticity and buoyant production respectively. Similarly (3.7) (or 3.9) contains the local rate of change of temperature, advection of temperature by the velocity field, and a thermal diffusion term. 3.2.1 Nondimensional Form of the Governing Equations Since the set of nondimensional parameters including length, velocity and time scales is not unique, there exists multiple possible representations of the equations in nondimensional form. Careful consideration of the physical problem under investigation may lead to clues as to an appropriate choice of scaling parameters. Typical choices for the length, velocity, and time scales used in the investigation of R-B convection, are given by ~ H 2 Length: [H], Velocity: v H , Time: (3.10) which is based on a viscous dissipative time scale and results in the introduction of the Prandtl number, and a Rayleigh or Grashof number. Appendix A discusses another possible scaling based on the atmospheric parameters of surface heat flux and convective velocity scale (w*). An estimate of the velocity scales based on (w*) is approximately an order of magnitude greater than a velocity scale based on (3.10). Since Chapter 3: Modified Rayleigh-Benard Convection ... 56 the selection of scaling parameters is often motivated by numerical limitations, it should be noted that computations based on the convective velocity scale (w*), will attain smaller values of the flow variables, but will necessarily be run for longer values of the nondimensional time when compared to a scaling based on the viscous time scale. However since the velocity scales differ only by a factor of 10, there is no computational reason to chose one form of the velocity scaling over the other. Thus the set of nondimensional parameters used in the current study are given by (3.10) and the following nondimensional variables are introduced H H ( T 0 - T O H 2 v The resulting nondimensionalized form of the governing equations for the streamfunction and temperature are d V 2 v F , 90 J 0 F , V 2 ¥ ) + G — - V 4 ¥ = 0 (3.11) dx 9X V 2 0 = O (3.12) ^ - J O F , 0 ) -g H where G - v 2 T 0 - T , V T o + T , y is a Grashof number and a = (v / K) is the Prandtl number . Thus the form of the solution obtained depends on the nondimensional parameters G and a . The numerical values of these parameters may be varied in order to identify criterion influencing the transition between various flow regimes. These transitions may be abrupt (e.g. initiation of convection) or may occur gradually over a range of the parameters (such as the transition from time-dependent to fully developed turbulent flow). See table 2.1 for convective regimes based on experimental observation. Chapter 3: Modified Rayleigh-Benard Convection . 3.3 The Governing Equations in (^,T() Space 57 3.3.1 Mapping from Physical Space to (£,,r\) Space Coordinate transformations are often introduced in order to simplify the application of boundary conditions. However, if the mapping is not orthogonal (as is the case for the current study), complications are moved from the boundaries to the governing equations in the sense of significant increases in the number of terms involved in such operators as the Laplacian. The modified R-B domain is mapped to the classical domain consisting of two flat plates, such that the boundaries defined in physical space (X,Z) by Z = 0 and Z = (1 + ecos(bX)), for all X , become the bounding planes r\ = 0 and r\ = 1, for all \\ in (^,r|) space (figure 3.3). zT — ' n t > X > £ Figure 3.3: Mapping of (X,Z) to (cj,T)) space. The mapping used is given by ^ = X ' r> = 71 <3"1 3) (l + ecos(bX)) where e is the amplitude of the forcing, b is the imposed wavenumber and the inverse mapping is naturally given by Chapter 3: Modified Rayleigh-Benard Convection . 58 Z = Ti(l + ecos(b£)). Based on the mapping given by (3.13) equations (3.11) and (3.12) become Temperature d® 1 j ( ¥ , 0 ) _Iv 2 0 = 0 3x (1 + ECOS(D£)) o (3.14) where V 2 0 = a 2 © a 2 © (1 l + r i 2 8 2 b 2 s in 2 (b^) aV dx\2\ (l + ecos(b^))2 + 2/2V f a 2 G +• 30 Tieb 2 •nebsin(bcj) (l + ecos(b£))) dr\ I (l + ecos(b$)) (2e + cos(bcj) - ecos 2 (b£)) Streamfunction 3 V 2 ¥ dx (1 + e cos(bcj)) J0F,V 2 v F) + G riebsin(bci) 3£ ' dr\\ (l + £cos(b^)) 30 30 + A' - V 4 v F = 0 (3.15) where V 2 v F is defined above and V 4 ¥ = 34*F * 1 4 1 (l + ecos(b4))4 +8 eb sin(bcj) (l + ecos(b^))3 | 3 4 xF | / ] d4x¥ ( riebsin(b^) ' + 3£ 4 + 3£ 3 3r|l (l + ecos(b£)) + 6-3 4 vF ( riebsin(b^) ^ 2 3^23ri21 (l + ecos(b^)) 44 d4y¥ r\eb sin(b^) 3£3ri31 (l + ecos(b^)) + -3V r(8bsin(bcj) (l + Ecos(b£)) Chapter 3: Modified Rayleigh-Benard Convection . 2 m / +4 a2 ^ •neb3 sin(b^) dSdnl (l + ecos(b^))3 (6e2 + 4ecos(bi;) - e 2 cos2(bi;) -1) +2 d3V an3 T|£b2 (l + Ecos(bc;))4 [4e sin 2 (b£) + 2e + cos(b^) - e cos2 (b£)] +4 ( „ u 2 „;„rut\ \ sb 2 sin(bc;) dri 2 I (l + £cos(bi;)) 4 [E sin(b^) + 2£ + cos(b^) - cos2 (bi;)] +(2£ + cos(bi;) - £ cos2 (bi;)) r|£b 2 (l + £cos(bi;)) a 3 y 3 \T/ / + 12 33*F T(£bsin(b^) ^( l + EcosCbi;)) + 6 a3 ^ dn3 ri£bsin(b?;) (l + £cos(b^)) +2-d4W 1 (l + £cos(b^)) 2 d4y¥ „ 3 4 v F ' + 2 r|£bsin(bc;) a 2 ^ ' _ a^an 3 U 1 + e c o s ( b ^)) + a 4rl ri£bsin(bc;) (l + £cos(bi;)) \ 2 4 \ d2¥( r | 2 e 2 b dri 2 I (l + £cos(bi;)) 4 [4 sin 2 (b^)(6£ 2 + 4£ cos(b^) - £ 2 cos2 (b£) -1) + 3(2£ + cos(b^) - £ cos2 a*F ( r|eb 4 A an I (l + £cos(bi;)) 5 ( E 3 COS 4(b^) - 1 1 £ 2 cos 3(b£) + (11£ - 20E 3)cos 2(bi;) - cos(bi;) + 8E(3E2 Chapter 3: Modified Rayleigh-Benard Convection ... 60 (Note that for 8 = 0 the (3.13) reduces to r\ - Z where 0 < r\ < 1 and 0 < Z < 1. However, for b = 0 but 8 * 0 , the (3.13) reduces to r| = Z / ( l + e) where 0 < r\ < 1 and 0 < Z < ( l + e). Therefore in the limit as b -» 0 (3.14) and (3.15), reduce to (3.11) and (3.12) with the addition of (1 + e) as a vertical scale factor.) From the definition of the streamfunction given by (3.5) the velocity components in mapped space become i (l + e cos(b£;) W(£,T|) = -W a^f er|bsin(bc;) A d t + d r i (l + ecos(bc;) (3.16) 3.3.2 Analytical Techniques In order to solve (3.14) and (3.15) two analytical techniques are used: • a weakly nonlinear perturbation technique which expands the variables in powers of e, • and an investigation into the fully nonlinear behavior of the equations resulting from the truncation of (3.14) and (3.15) to first order in e. A weakly nonlinear perturbation analysis is used to investigate the behavior of the solutions resulting from an expansion of the streamfunction and temperature to first order in e. The disadvantages of a weakly nonlinear analysis were discussed in section 2.1. Recall that the regime of applicability of the solutions obtained using a weakly nonlinear analysis is limited to Chapter 3: Modified Rayleigh-Benard Convection ... 61 near critical values of the Grashof number. Therefore it is anticipated that a weakly nonlinear perturbation analysis will provide useful information with respect to the early stages of flow development. Explicitly, *F and 0 are assumed to be the form Y = y 0 (x ,^Ti) + e Y I ( T , t T l ) + . . . 0 = 0 O (x, TI) + e©, (x, ri)+... (3.17) where the equations governing the leading terms are denoted zeroth order equations (with subscript '0'), and the equations governing the second term in the expansions are denoted first order equations (with subscript T ) . The expansion given by (3.17) differs from 'standard' perturbation analysis techniques used to study the classical R-B convection problem (e.g. Busse, 1967) since the standard approach expands the solution about the conductive base state, and the first order terms represent the secondary roll circulation (see sections 2.1 and 2.2). In contrast, here the zeroth order equations represent the fully nonlinear classical R-B convection problem and naturally contain both the conductive and convective solutions depending on the value of the Rayleigh (Grashof) number. Unlike the standard approach, the current study does not impose a small amplitude restriction on the roll circulation and is fully nonlinear. In the current study the first order terms represent the circulation which is effected by the warping of the upper surface. The expansion (3.17) does however, assume that the first order solution is small compared to the zeroth order solution, and that second order effects of the corrugated surface are negligible. Although the zeroth order equations are nonlinear (section 3.3.3), the first order equations are linear (section 3.3.4) thus suggesting that the perturbation technique may not model Chapter 3: Modified Rayleigh-Benard Convection ... 62 accurately the nonlinear flow dynamics. This lead to the consideration of a simplified version of (3.14) and (3.15) where the expansion of the mapping is truncated to first order in the amplitude of the warping. This truncated system of equations was felt to retain more of the crucial physics of the original system (section 3.3.5). The limitations of the perturbation analysis and comparison of results obtained from both analytical techniques will be discussed and presented in more detail in relevant sections. 3.3.3 Zeroth Order Equations Upon mapping of the coordinates from (X,Z) to (£,T|) space, and substitution of the expansions for the streamfunction and temperature given by (3.17), the equations governing the zeroth order variables lead to the following expression for the temperature ^ - J O F o , 0 o ) - i v 2 0 o = O ox a (3.18) and streamfunction 8V 2 ¥ , dx 0 - J (*F 0 ,V 2 ^ 0 ) + G ^ - V > n =0 (3.19) where V 2 = a2 a2 a^2 an2 and V 4 = a4 „ a4 + 2 + -a^ 4 aifan2 an Note that the zeroth order equations are both nonlinear and homogeneous. Comparing (3.18) and (3.19) with (3.11) and (3.12), it can be seen that the effects of forcing at the upper surface does not appear in the equations governing the behavior of the zeroth order variables. Thus (3.18) and (3.19) correspond to the governing equations for the nonlinear development of a Chapter 3: Modified Rayleigh-Benard Convection ... 63 fluid confined between two flat, rigid plates of infinite horizontal extent, i.e. classical R-B convection. 3.3.4 First Order Equations The equations governing the first order variables, are the following: for temperature, ^ - - J ^ ^ G o ) - ^ , © , ) - - ^ 2 © , = H O F O , 0 O ) ax a (3.20) with H ( v F o , 0 o ) = - V 0 o - c o s ( b ^ ) J ( v F o , 0 o ) and d2 d2 d 2r|b sin(b^) — — 2 cos(b^) —j + nb 2 cos(b^)— dqdri dr\ dn. and for the streamfunction, , dx , - J ( Y 1 , V 2 Y o ) - J ( T o , V 2 Y 1 ) - V 4 Y 1 + G ^ - = Q ( T o , 0 o ) (3.21) with Q(vp o ,0 o ) = - ^ _ A J +J (Y 0 ,V7Y 0 ) -cos (b^)J (Y 0 ,V 2 Y 0 ) +4 rubsin(bcj) V 2 T n - cos(bcj) V 2 Y n 6 T 1 b 2 c o s ( b 0 ^ g - + 8 b s i n ( b ^ ) ^ - + 2 T 1 b 2 c o s ( b ^ ^ dc, dr\ d^dr] dr\ 4nb 3 s i n ( b ^ ) ^ - - 4 b 2 c o s ( b ^ ) | ^ nb 4 c o s ( b ^ ) ^ - G n b s i n ( b ^ ) ^ dri dri Chapter 3: Modified Rayleigh-Benard Convection ... 64 Examining the terms on the LHS of (3.20) and (3.21), we note that the first order equations are linear in the first order variables. Since the solution of the zeroth order equations for a given Prandtl and Grashof number are known, the terms associated with the Jacobian do not result in nonlinear terms, but in terms involving nonconstant coefficients. Also note that although the forcing functions H O F O , 0 O ) and QOF0,eo), which arise on the RHS of (3.20) and (3.21) are complicated and nonlinear, they depend only on the zeroth order variables. Therefore, in order to solve the first order equations, the solution of the *F0 and 0 O satisfying (3.18) and (3.19) must be known. 3.3.5 Truncated Equations As noted in sections 3.3.3 and 3.3.4, the perturbation analysis results in nonlinear, homogeneous equations for the zeroth order terms, and linear, forced equations for the first order variables. Thus, the form of the perturbation expansion chosen (3.17) results in the implicit assumption that the solution is unaffected by the forcing at the upper surface to zeroth order in e, and therefore does not allow for modification of the basic (zeroth order) solution in response to the upper forcing. Also, the linear response of the first order terms, which are forced by an assumed known base state (i.e. solution of the zeroth order equations) does not allow for the transfer of energy between modes of the first order equations or between first order modes and the zeroth order base state. Therefore, in addition to the investigation into the response of the zeroth and first order equations, we will also consider the effects of the forcing on a truncated version of the complete Chapter 3: Modified Rayleigh-Benard Convection ... 65 set of governing equations (to first order in epsilon). The truncated equations are the following the temperature ^®- = J eF,0) + - V 2 0 + eHOF,0) (3.22) dX O and streamfunction d_ dx de (V'+eV)*? = J0F ,V2 v F) + V 4 v F + G — + eQOP,0) (3.23) J ^ where fR^F,©) and Q(*F,0) are as defined below equations (3.20) and (3.21) respectively. Comparing (3.22) and (3.23) with equations (3.18) and (3.19) which govern the behavior of the classical R-B problem, the truncated system of equations incorporates the additional terms H( v F,0) and Q( l F,0) on the RHS of the temperature and streamfunction equations respectively, where the amplitude of the terms is regulated by the magnitude of e. Thus, the effect of the mapping of the governing system of equations is the addition of terms which are complicated expressions in *F and 0. 3.4 Boundary Conditions For all of the equations under consideration it was assumed that the surfaces at r| = 0 and T| = 1, were rigid, of a constant temperature, and that no stress is exerted by the flow at the boundaries (i.e. free-slip flow conditions). Explicitly the boundary conditions for the zeroth order variables, and the truncated equations are given by dW d2W dy d 2 V F r l = 0 : ^ = 0'4p^ = 0 a n d 0 o = 1 TI = 1: ^ = O , ^ = 0 and e 0 = - l . Chapter 3: Modified Rayleigh-Benard Convection 66 and the following for the first order variables T| = 0: = 0 d2% = 0 and 0 , = 0 d2% TI = 1: d?, = 0, an2 = 0 and 0 , = 0 . 3.5 Fourier Integral Analysis As noted in section 2.5.2, in order to avoid the restriction of prescribing periodicity over any predetermined length scale, a Fourier integral representation for the streamfunction and temperature in the horizontal is used following Getling (1983). The vertical flow dependence is assumed to be adequately represented by a finite Fourier series. The transformation of the equations was done via a two step process. The first step involves eliminating the vertical dependence on n by expanding the flow variables (represented by F(x, cj, T|)) in terms of a Fourier sine series: i.e. where Fk(t,cj) = 21 F(t, cj, r\) sin(k7tn)dri are the Fourier sine coefficients. The expansions for the flow variables are then substituted into the terms in the governing equation. Using Galerkin's method, the orthogonality of the sinusoidal basis functions is exploited in order to obtain an equation governing the temporal development of the Fourier sine coefficients. The second step is to transform the resulting equations to wavenumber space (a) using a Fourier integral transform (3.24) k=l 0 Chapter 3: Modified Rayleigh-Benard Convection ... 67 fk(x,a) = JPfcCc^e-^dS (3.25) where fk(x,a) is the Fourier transform of F k(x,cj). The inverse transform is given by F k ( x ^ ) = ^ j f k (x ,a )e^da . (3.26) Note that the Fourier coefficients fk(x,a) are written f a k for brevity. Requiring F(x,c;,r|) to be real, the solution in (£,r)) space may be obtained using (3.26) with F k ( x , = J - J [faRk cos(a$) - fa',k an(a$)]da (3.27) 27U_ where the superscript 'R ' and T denote real and imaginary component respectively. Also f i = f - . . k a n d f. ' k=-f-.. k- ( 3 - 2 8 ) 3.5.1 Zeroth Order Equations in Wavenumber Space The form of the expansion for and 0 O which satisfy the boundary conditions given explicitly in section 3.4.1, may be written V0M,r\) = £ ¥ k 0 ( x ^ ) s i n ( k 7 r n ) (3.29) k=l -I'M 0o(x,c;,Ti) = (l-2Ti) + X© o k(i:^)sin(k7rri) (3.30) k=l where from (3.26) Y K ° ( T , ^ ) = — ? V a k e ^ d a and 0 ° ( T , ^ ) = — ?0° k e i a ^da . Note that the 2K ^ ' 2n ^ ' constant temperature boundary conditions are satisfied with no restrictions on 0 k . Chapter 3: Modified Rayleigh-Benard Convection ... 68 Upon substitution of the above expansion for *F0 and 0 O given by (3.27) and (3.28) into the governing equations (3.18) and (3.19), and converting to wavenumber space, results in the following integrodifferential equations for 0 a k and \ | / a k respectively; Temperature 3ft" 1 ° ° ° ° ° ° ^ = -^Ylke:,k-i2a<k+lXIJ[cU,m,k]¥L,,me^nda' (3.31) OX (J 71 M = 1 N = I _^ Streamfunction - Y * k ^ = - i a G 9 ° k + Y:,k<k - i X I jYj.B[cU. n, k]v?.-..xm< nda' (3.32) ° " n m=l n=l T E R M (1) withy 2 , y 4 k and C a a . n m k defined in appendix B by (B.2). Note that the subscripts of the coefficients signify their dependence on vertical mode numbers (n,m or k) and/or horizontal wavenumber (a, or a'). The Fourier analysis has resulted in the transformation of the nonlinear terms from differential form to integral form (TERM (1)) where the integration is over the entire wavenumber spectrum. The double summation arises from the product of the vertical mode dependence of the Fourier coefficients, where the orthogonality of the basis functions are accounted for within the coefficient C a a . n m k . The equations governing the real and imaginary components of the streamfunction and temperature are given explicitly in appendix B. Chapter 3: Modified Rayleigh-Benard Convection ... 69 3.5.2 First Order Equations in Wavenumber Space The expansions for x¥ 1 and 0 j which satisfy the boundary conditions of section 3.4.1 are given by *F, (T , Tl) = X (x, $) sin(]acn) (3.33) k=l 0 1 (T,^,Ti) = X© k (x , c> in (k7 tn ) (3.34) k=l where Y ' ( x ^ ) = j - } < k e ^ d a and ©[(x,^) = j - j e a k e ^ d a . For the analysis of the first order equations, it is assumed that the zeroth order solution is known and is represented by the solution to the equations presented in section 3.5.1, and therefore the forcing (which is only a function of the zeroth order solution) may be evaluated explicitly. Recalling the governing equations for 0 , and  XVX given by (3.20) and (3.21), substituting in the above expansions, and transforming to wavenumber space, we obtain the following governing equations for 6 a k and \|/ a k Temperature ^ - = - 1 Y l k e : , k - i 2 a ¥ I , k dx a JcI,,,„,m,k[vU),me°,„ + < ^ n e , , , ] d a ' + 2 h « k , e ° k ) (3.35) m=l n=l . where Chapter 3: Modified Rayleigh-Benard Convection . 70 2 i i L [ § ( a - b ) + 5(a + b)] + i ( k7 t ) 2 [ 0 ^ b ) , k + 0 | , a + b ) , k ] T E R M (1) T E R M (2) +"I(n7C)Pn,k[(a" b/2)0°_ b ),„ - (a + b/2)0° a + b ) , n ] O n=1 A J[^a,a',n,m,kV(a-a'-b),m 4 '^a,a',n,m,kV(a-a'+b),m^a , ,n^ a m =l n=l _ „ T E R M (3) Streamfunction ,2 dt|/a_k _ . _ ^ a 1 _ _ 4 a,k ax - iaG0: , k + Y : , k xi / a,k jYa ,„[cI , , , n ,m ,k][¥U),m¥°,„ + <-,),,„¥:,•,,>• + 2 q « k , 0 ° k ) (3.36) ft m=l n=l where q ( ¥ : , k , e : , k ) = iGb 2 L 2 L [8(a + b) - 8(a - b)] + £(n7t) |3 n , k I [eJ._b),n - 0° + b ) , n ] n=l T E R M (4) - ( k 7 t ) 2 Y 2 i k [ < a _ b ) , k + < +b),k] +X^ n 7 t ' ) P n , k [ a 2 ¥ | ) a - b ) , n + a+2 ¥(l+b),n ] n=l + I K V(a-a'-b),m + « 3 V(a-a' +b),m ]Va ' , n da' m=l n-1 . with p „ j k , Xk, n .m , X k , „ , m > C 2 ; . n , m , k , C 2 ; . n i m , k > a " , a*, a" and a* defined by (B.4). Although the expressions for h(\ | /" k ,0 a k ) and q ( t | / a k , 0 a k ) are complicated, some insight into the behavior of the solution of the first order variables may be obtained by considering Chapter 3: Modified Rayleigh-Benard Convection ... 71 TERMS (1), (2) and (3) noted as part of the expression for h ( \ | / ° k , e ° k ) , and T E R M (4) noted as part of the expression for q C V a . k ' ^ a . k ) • The delta function of T E R M (1) arises as a result of the product of the vertical component of the velocity field and the linear portion of the temperature expansion for the zeroth order variables (i.e. the (1 - 2r|) term). This term is analogous to the Fourier transform equivalent of the expression w' (d© / dz) (where w' is the vertical component of the perturbation velocity, and (d© / dz) is the mean Vertical potential temperature gradient), which arises from the advection term in the thermodynamic equation, following a (typical) Reynolds decomposition of the flow variables. Since T E R M (1) is real, it will form part of the forcing function for the equation governing the real component of the temperature only, where we note that the value of the delta function has a zero value for all wavenumbers other than (a = b) or (a = - b ) . Since T E R M (1) has a nonzero value for all b with a magnitude greater than zero, this implies that there exists a forcing at these values of the wavenumber spectrum, whose magnitude is dependent on the square of the imposed wavenumber 'b'. Similar to T E R M (1), the delta function associated with T E R M (4) has nonzero values only for wavenumbers identical to the imposed wavenumber 'b'. However, this term arises as a result of the coordinate transformation of the buoyancy term G ( 3 0 / d X ) of (3.11) to mapped space. Since T E R M (4) is imaginary, upon separation of the equations into real and imaginary components, the delta function will form part of the forcing function for the imaginary component only. Of particular interest is the fact that this term has a finite and nonzero value for all (bounded) values of G and Ibl greater than zero. (See appendix B for the equations governing the real and imaginary components of the streamfunction and temperature.) Chapter 3: Modified Rayleigh-Benard Convection ... 72 Terms of form (3) are nonlinear and arise in association with the mapping of the advection terms. Finally, T E R M (2) is typical of the remaining terms of h ( \ | / ° k , e ° k ) (or q ( \ | / ° k , 8 ° k ) ) where the dependence on the solution of the zeroth order equations is shifted both left and right in wavenumber space. These terms have a zero value for subcritical values of the Grashof number (i.e. within the conductive regime) and a nonzero value for supercritical values of the Grashof number (i.e. within the convective regime). Section 3.6.2 will present a linear analysis of the equations designed to investigate in detail the response of the first order equations to a linearized forcing. 3.5.3 Truncated Equations in Wavenumber Space Consider the LHS of (3.23) given by — ( V 2 +eY)x¥. Since we are looking for the steady-state dx 3V 2 ¥ solution and the parameter e is assumed to be small, this term is approximated as . dx Justification for such a simplification is outlined in appendix C through the use of scaling arguments and post-consideration of numerical results. However it should be noted that the approximation was motivated by the resulting simplification in the numerical analysis. If the e term is included, then once converted to wavenumber space the LHS would result in a system of coupled equations which would need to be solved simultaneously rather than equations governing the Fourier coefficients which may be solved using explicit methods. Chapter 3: Modified Rayleigh-Benard Convection ... 73 The expansion used for the flow variables associated with the truncated equations, is the same as that used for the zeroth order variables, i.e. (3.29) and (3.30) Transforming (3.22) and (3.23) to wavenumber space leads to the following Temperature 30 1 • ^ + - Y a ^ e a , k + i 2 a ¥ a , k - i X X j ^ (3.37) CT CJ 71 M = 1 N = 1 _^ Streamfunction ~ y ; U ^ ^ + i a G e a , k - Y : , k ¥ a , k + ^ =82q(x(/ a k ,0 a k ) °t >l m=l n=l _«o (3.38) whereh(\ | / a k ,0 a k ) and q (y a i k , 0 a > k ) > are defined by (3.35) and (3.36). The coefficients are defined by (B.2) and (B.4). 3.6 Analytical Solutions The following sections present the results of an analytical investigation into the behavior of a linearized version of the zeroth order, first order, and truncated equations in order to provide insight into the anticipated response of the fully nonlinear system. 3.6.1 Zeroth Order Equations As discussed in section 2.2 the solution to the linearized version of the zeroth order equations was first presented by Rayleigh (1918). As results are well documented, they will be summarized here Chapter 3: Modified Rayleigh-Benard Convection ... 74 by noting the relationship between the Grashof number and the wavenumber of the only nonzero solution which is given by » = f ^ . (3.39) 2a a The value of the wavenumber corresponding to the critical Grashof number (G c ) may be found by minimizing the expression for G with respect to the wavenumber. This results in a minimum wavenumber a c - 2.221 and G c = [(3/ 2 ) 3 7 t 4 a _ 1 ] = 328.75a"'. 3.6.2 First Order Equations with a Linearized Forcing In order to investigate the behavior of the first order equations to a linearized version of the forcing, it is assumed that • the zeroth order solution for Y 0 and 0 O is known and of the form Y 0 = V ( a „ ) sin(TCTi) 0 O = (1 - 2TI) + 0(a o) sin(7m.) (3.40) where 6(a 0) = A R + i A 1 , V|/(a0) = B R + i B 1 and the real and imaginary coefficients are related by A R = [2aa/(a 2 + 7 t 2 ) ] B I and A 1 = [-2aa/(a 2 +TT.2)]BR ; • the Grashof number and preferred wavenumber (a 0) of the zeroth order solution are related by (3.39); and • the first order variables contain only the first vertical mode. (Note that the assumed form for the solution of the zeroth order equations is consistent with a linear solution containing one vertical mode, and is motivated by analytical simplicity. Such a solution will be most valid for values of the Grashof number near critical.) Chapter 3: Modified Rayleigh-Benard Convection . 75 Introducing the vector of unknowns y = [flf , 0 a , \ | / R , the linearized form of the steady-state, first order equations may be written in matrix notation as Ay = F (3.41) with A = ( a 2 + 7 T 2 ) 0 0 0 ( a 2 + 7 X 2 ) 2ao 0 -aG - ( a 2 + 7 t 2 ) 2 2ac> 0 0 aG 0 0 2 x 2 -(a z + 7 X Z ) and F = f, + f2 + f3 + f4 (3.42) The components of the forcing vector are defined such that and - 4 b 2 0 0 4Gb if a = b else f. = 0 , -4b 2 ' 0 0 -4Gb if a = -b else f2 = 0 2 A R 2A 1 n2 — ( a D + b /2 ) 4 ° 7 X 2 - ^ ( a 0 + b / 2 ) 2 A 1 — - 2 B R - 2 A R ^ - 2 B I ' (a„+b) ^ 2l(a 0+b) - 2 . , 2 1 ft Y r „ +M + — '(a„+b) ^ 2 l(a0+b) if (a - b) = a 0 else f3 = 0 2A 2A 1 . 2 b Kl+-(a - b / 2 ) 4 0 . 2 b 7 t / + - ( a 0 - b / 2 ) _ 2 A I ^ - 2 B R 8 2 A R ^ - 2 B ' 8 7 1 Y(..-b, + 4 « 2 | ( 1 . _ b ) 1 . 2 . . 2 r Y f.„- b)+-a 4 2 l (a 0 -b) if (a + b) = a 0 else f4 = 6 Chapter 3: Modified Rayleigh-Benard Convection ... 76 - | b 4 -2ba 0 (a 0 2 + 7 t 2 ) - 2 b 3 a 0 - b V - 3 b 2 a 2 - - b 4 +2ba b(a 2 +7t2) + 2b 3 a 0 - b V - 3 b 2 a 2 In order to simplify analysis of the first order equations, the following two cases were considered; the conductive zeroth order state with G < G c ; and the convective zeroth order state with G > G c . Use of the superposition principle (for both of these cases) allows for the consideration of each of the forcing components individually, in order to develop the net, first order response. The Conductive Zeroth Order State Consider G < G c and therefore the zeroth order solution is Y 0 = 0 , 0 O = ( l - 2 r | ) , and (3.41) reduces to Ay = f, + f 2. For a ^ ± b , f, = f2 = 0 and (3.41) reduces to the linear zeroth order equations consider in section 3.6.1. Since G < G c this suggests that all Fourier coefficients for wavenumbers other than those corresponding to a = ±b are zero. For a = b (3.41) reduces to A | a = b y = f, (3.43) where the subscript denotes that the coefficient matrix A is evaluated at a = b , and A and f, are defined by (3.42). The solution vector may then be found using y = A ' ^ f , where A J n v is the where: a , 2 (a„+b) a l(a„-b) Chapter 3: Modified Rayleigh-Benard Convection ... 77 inverse of A | a = b . However care must be taken to ensure that the coefficient matrix is invertable, and since the governing equations may be written in coupled pairs for y i =[9 b R , i | / b ] and y 2 = [ e b , t ) / b R ] (3.44) (3.43) may also be written A 2 3 y . = [0>0]T> and A 1 4 y 2 = 4b[-b,G] T where ^23 _ ( b 2 + 7 t 2 ) 2ba - b G - ( b 2 + 7 t 2 ) 2 and A 14 ( b 2 + 7 t 2 ) -2ba bG - ( b 2 + 7 c 2 ) 2 Solving for yx with G < G c gives y, = [0,0]T . Similarly, the solution for y 2 may be found by inverting the coefficient matrix A i 4 giving y 2 - x b [ ( b 2 + 7 t 2 ) 2 + 2 a G ] G(2b 2 + 7 t 2 ) (3.45) with X = [(b2 + 7 t 2 ) 3 - 2b 2 aG]. If b satisfies (3.39) then X = 0, and forcing vector results in an unbounded amplitude for the real component of temperature and the imaginary component of the streamfunction. This is due to the system being forced at its preferred ('resonant') wavenumber. Also for b ^ 0 then y 2 * [0,0]T. Therefore, for all Grashof numbers greater than zero, a first order circulation is induced for all imposed wavenumbers. The effect of the horizontally varying temperature field resulting from the upper corrugated surface may be viewed by analogy to the problem of flow between vertical plates of uniform, but different, temperatures. Appendix D elaborates on this comparison between a corrugated surface and flow between vertical bounding Chapter 3: Modified Rayleigh-Benard Convection ... 78 planes and presents a proof that the horizontal temperature gradient results in initiation of a circulation. Finally, it can be shown that for a = -b and using (3.28), the same solution as for a = b is obtained. The Convective Zeroth Order State For G > G c with the zeroth order solution approximated as outlined in the assumptions, a number of subcases are considered: For a ^ b , a # - b , (a - b) ^ a 0 , and (a + b) ^ a 0 , (3.41) reduces to the homogenous form for which the only nontrivial solution is for the wavenumber satisfying (3.39). The solution for a = ±b has been considered and the solution is given by (3.45) and yl =[0,0] where y, and y 2 are defined by (3.44). For (a -b) = a 0 , (3.41) becomes A| a = a + b y = f3 where A and f3 are defined by (3.42). The solution vector may then be found by inverting the coefficient matrix such that y(a 0 + b) = A f f3 where A 3 n v is the inverse of A | a = a o + b . Similarly, for (a + b) = a 0 , (3.41) becomes A| a = a _ by = f4 where A and f4 are as defined by (3.42). As for the case that (a - b) = a 0 , the solution vector may be found by inverting the coefficient matrix, y(a 0 - b) = A 4 v f 4 where A 4 n v is the inverse of A| a = a _ b. Therefore for given values of G, b, a,,, and amplitudes of the zeroth order solution, we may invert the coefficient matrices A3 and A 4, and multiplying by the forcing vector, obtain the magnitude of the first order coefficients for the temperature and streamfunction. Chapter 3: Modified Rayleigh-Benard Convection . 79 Summary Linear analysis of the first order equations has provided interesting insights into the anticipated response of the physical system. For subcritical values of the Grashof number the only nontrivial solution is given by (3.45) for a = ± b , where b is the imposed wavenumber. Since a nonzero solution occurs for all b ^  0 with 8 ^ 0 and G > 0, there is no stable base state for the modified R-B convection model in contrast to the classical problem. For supercritical values of the Grashof number associated with a convective zeroth order state, the solution vector governing the first variables is given by " e R ~ ° ( a 0 + b ) ° ( a 0 - b ) e i - < + + al °(a„+b) ¥ ( a 0 + b ) + ft1 U ( a „ - b ) V d - b ) _ ¥ ( a „ + b ) _ ¥ ( a 0 - b ) . where a 0 satisfies (3.38). Also, "6?" et 4b ¥ b X b(b 2 + 7 t 2 ) 2 +2aG' 0 0 G ( 2 b 2 + 7 t 2 ) and o R ° ( a 0 + b ) el. (a0+b) R (a0+b) ¥ ( a 0 + b ) . A^ n v f 3 with o R D ( a „ - b ) ° ( a „ - b ) ¥ ( l - b ) - A i n vf - A 4 i 4 . (Note that ^ is the wavenumber corresponding to the assumed zeroth order solution, i.e. The first order equations contains the same solution as the zeroth order equations, and would be the only solution if b = 0.) Chapter 3: Modified Rayleigh-Benard Convection ... 80 Note that as b -> a 0 and G -» G c the magnitude of the streamfunction increases. This is not surprising, since forcing at the preferred scale will result in resonance and subsequent unbounded growth under a linear framework. (See appendix E for analytically predicted values of the streamfunction and temperature components.) 3.6.3 Linearized Form of the Truncated Equations In order to investigate the behavior of the linearized, steady-state, truncated equations, it is assumed that that the solution consists of a single mode of the form i i ) = x¥l sin(TXTi) 0 & TI) = (1 - 2TI) + 0 , sin(Trn) where the subscript T denotes the first vertical mode. (The subscript '1 ' is subsequently dropped for brevity.) Separating into real and imaginary components (3.37) and (3.38) become Temperature y2eaR - 2 a a ¥ a - e [ A : G R _ b ) + A+eR+b)] = -4b 2 e[8(a-b) + 5(a + b)] y20a + 2aa ¥ a R - e[ArG| a_ b ) + A ^ U ] = 0 (3.46) Streamfunction - y V f -aG0l +e [B: ¥ R _ b ) + B : ¥ R + b ) ] - E y [ ( e | 8 . b ) -ej 1 + b ) )] = 0 (3.47) - Y > | + aG9 R + e [ B 7 V U + B [ v U ] + e ^ [ ( 0 U -eR+b))] = 4Gbe[8(a-b)-5(a + b)] h b 1 where A7=2 (7t 2 ( a - b / 2 ) ) , A * = 2(TX2 + - ( a + b /2) ) , B," = 2(7t 2(a 2 +7X 2) +—a 2) and 4 4 4 B^ = 2 ( 7 X 2 ( a 2 + 7 X 2 ) + | a 2 ) . Chapter 3: Modified Rayleigh-Benard Convection . 81 Note that (3.46) and (3.47) may be written as pairs of coupled equations in (6 a,\)/R) and (0f,\j/a). If it assumed that for all wavenumbers Ial> a m a x 0 ^ k = V^L.* = vfLx.k = 0 ' then (3.46) and (3.47) may be written in matrix form as Ay, =0 B y 2 = F (3.48) where Y?..„B). - E A * | ( . . N B ) — ...,2a(a-nb), 0,... - e A > l ( a - ( „ - 1 ) b ) ' YU-IW - e A : | ( a - ( „ - , ) b ) ' - ° ' 2a(a-(n-l)b), 0,... - - e A > l ( a - b ) ' - e A l | ( a - b ) ' ' 0 ' 2 ° ( a - b ) v 0,... • • • - e A r | w . Y(a)> - £ A i | ( a ) ' >0> 2 G a , ° — • • • ' - e A i | ( a + b ) ' Y U - - e A ' L ) ' ' ° ' 2o-(a + b), 0,... . . . - e A j i , ylUnUM, -eA+l ,. ,0, 2a(a + (n-l)b), 0 1 l(a+(n-l)b)' '(a+(n-l)b)' 1 l(a+(n-l)b) ' V V J " ~ e A i " l ( a + „ b ) ' T ( . + n b ) . - --.0, 2rj(a + nb) -G(a-nb), eGb/4, 0, -ylDb), e B * | ( a n b ) , 0,... -eGb /4 , -G(a- (n- l )b) , eGb/4, , eB: | , - Y * eB*| v v ' ' ' ' 1 l(a-(n-l)b)' '(a-nb)' 1 l(a-(n-l)b)' . . . - e G b / 4 , -G(a -b ) , eGb/4, » e Br| ( a. b )» -JU)> e B i " l ( . . b ) ' - - • . . . ,-eGb/4, - G a , eGb/4, > e B i | ( a ) ' -Y1>. e B^| ( a )>--eGb /4 , - G ( a + b), eGb/4, ~lUv eB*L+b)— . . . ,-eGb/4, - G ( a + (m-l)b), eGb/4, ,eB~l , -yi+lmUM, eB+l ' ' V V / / ' > > 1 l(a+(m-l)b)' ' (a+(m-l)b)' 1 l(a+(m-l)b) . . . - e G b / 4 , - G ( a + mb)... 0, e B : | ( a + m b ) , -y^+m Chapter 3: Modified Rayleigh-Benard Convection ... and it has been further assumed that I (a - (n + l)b)l> a m a x and I (a + (m + l)b)l> a m a x . Also B = "YU)» - e A i | ( a - „ b ) ' - ...,-2o(a-nb), 0,... - e A r | ^ . , w « y l ^ r - e A - L - W ' ' ° ' -2a(a-(n- l )b) , 0,... • • • - e A : | ( a ) , y 2 a ) , - e A ^ | ( a ) ' ' ° ' _ 2 f 7 a ' ° ' -• • • ' - e A i | ( a + b ) ' Y U ) . -e A . | ( a + b ) ' > ° > -2^ (a + b), 0,... ••••-E A»U1 W' YU .W - ^ L W ' ° ' 2a(a + (n-l)b), 0 _ e A i | ( a + „ b ) ' Y(a+nb)' ,0, 2a(a + nb) G(a-nb), eGb/4, 0, - y 4 a . n b ) , eB+| ( a n b ), 0,... eGb/4, G(a-(n-l)b), - e G b / 4 , . e B : ^ , - y 4 a . n b ) , e B ; ^ , . . . e G b / 4 , G(a-b), - e G b / 4 , >eBi"|(l_b)> -Y(a-„v £ B i | ( a . b ) > -. . . ,eGb/4, Ga, - e G b / 4 , ,eB:| ( a ), - y 4 a ) , eB+| (a)>... . . . ,eGb/4, G(a + b), - e G b / 4 , >eBi"|(a+b)> - Y j + b ) » e B i | ( a + b ) > -. . . ,eGb/4, G(a + (m-l)b), - e G b / 4 , . e B : ^ i ) f e ) , - y 4 a + ( m , ) b ) , ^l+(m_m . . . ,eGb/4, G(a + mb)... .... 0, e B : | ( a + m b ) , - y 4 a + m b ) Chapter 3: Modified Rayleigh-Benard Convection . 83 and '(a-nb) '(a-(n-l)b) '(a-b) (a) 8' (a+b) I (a+(m-l)b) I (a+mb) ¥ ( a - n b ) ¥ ( a - ( n - l ) b ) (a) < + b ) R (a+(m-l)b) V(a+mb) " A R ° ( a - n b ) ° ( a - ( n - l ) b ) 6<a -b) e(Ra) ° ( a + b ) 0( Ra +( m-l)b) A R °(a+mb) I (a-nb) I (a-(n-l)b) V!.-b) Via) ¥ | a + b ) I a+(m-l)b) I (a+mb) Finally F = f I + f 2 where the only nonzero elements of fj correspond to 9 b and \ j / b , and has values of -4eb 2 , and 4eGb respectively. Also the only nonzero elements of f2 correspond to 9 R b and \ | / l b and has values of -4eb 2 , and -4eGb respectively. Following the linear analysis of the first order equations presented in section 3.9.2, the nature of the induced solution for nonzero values of the Grashof number is investigated. A Chapter 3: Modified Rayleigh-Benard Convection ... 84 comparison of the response predicted by the truncated equations with those obtained from the perturbation analysis is presented. The Conductive Regime Consider G < G C . Since (3.46) and (3.47) are complicated, the linear analysis focuses on wavenumbers associated with the imposed wavenumber 'b' and its harmonics (as the harmonics are coupled through the 'shifted' terms of (3.46)). For a = b the governing equations for the real component of temperature and the imaginary component of the streamfunction may be written in matrix notation as [ B y 2 L b = f - - (3-49) Since0 R a = 0 R , -0 r_ a = 0 a , \ | / R a = \|/Rand - \ | / a = \|/ a the equations may be simplified to contain only the components associated with 0, b, 2b, 3b ...mb, where it is assumed that (m + l)b > a m a x . Using M A P L E , the solution of (3.49) for the special case of m = 3 for G = 1 and G = 250 with G = l , b = l , and e = 0.01, was found. Results indicate nonzero values of the streamfunction and temperature, not only at values equal to the imposed wavenumber, but at the harmonics of the wavenumber as well. This implies that the warping of the isotherms associated with the corrugated upper lid results in multiple-scaling and not just single-scales as predicted by the perturbation analysis of section 3.6.2 for subcritical values of the Grashof number. It should be noted however, that the nonzero values at the harmonics of 'b' would arise in association with higher order terms of the perturbation analysis, which have been neglected. Details of the numerically evaluated cases obtained from the analytical and numerical models are presented in appendix E. Chapter 3: Modified Rayleigh-Benard Convection ... 85 3.6.4 Discussion of the Analytical Results A limited linear analysis has been conducted on the first order and truncated sets of equations. The objective of the analysis was to determine the linear effect of the corrugated upper surface on flow dynamics. In particular, the analysis focused on scales which are introduced into flow dynamics as a result of the warping of the upper surface. Results from the first order equations indicate that a circulation is induced at the imposed wavenumber for all nonzero values of the Grashof number. For G < G c the induced circulation was found to be of a single scale. For G > G c the warping introduces scales at the imposed wavenumber b, (a 0 +b), and ( a 0 - b ) where a<j is the wavenumber associated with flow development for classical R-B convection (i.e. the zeroth order equations). Consideration of subcritical values of the Grashof number in association with the truncated equations resulted in an induced circulation not only at 'b', but at the harmonics of 'b' as well. Thus, in contrast to the first order results the truncated equations predict multiple scaling for G < G C . The discrepancy between results from the first order and truncated equations with respect to the scales which are made available for nonlinear interaction with the flow will be shown to be significant (section 3.8). Furthermore, interpretation of results (theoretical and numerical) based on the perturbation analysis will be shown to be misleading when compared with results from the truncated equations. Chapter 3: Modified Rayleigh-Benard Convection ... 86 3.7 Numerical Analysis and Implementation This sections outlines how solutions of the integrodifferential equations were obtained numerically using standard numerical techniques (Press et ah, 1992). 3.7.1 Transformation of the Analytical Equations to Numerical Representation In order to obtain an approximate, numerical solution of the zeroth order ((3.31) and (3.32)), first order ((3.35) and (3.36)) and truncated ((3.37) and (3.38)) set of equations, the analytical system of governing equations for the real and imaginary components of the streamfunction and temperature (which assume the summation over all vertical modes and integration over the entire wavenumber spectrum) must be truncated at some finite value. Thus the expansions for the streamfunction and temperature are truncated such that and it is assumed that N and A are chosen large enough to ensure that the flow has as much freedom as numerically possible, thereby allowing for the 'natural' selection of the preferred wavenumber. Since numerical integration requires the use of discrete points, the analytical continuum in wavenumber space must also be discretized. Although this forces the solution to attain a nodal value, the resolution is only restricted by the number of nodes chosen to represent the wavenumber spectrum. The number of nodes may be chosen in order to obtain as high a degree of resolution as desired, at the expense of computational speed. Chapter 3: Modified Rayleigh-Benard Convection ... 87 The relationships between the real and imaginary components for positive and negative values of the wavenumber given by (3.28) allows for consideration of only positive wavenumbers. Example Consider the nonlinear term associated with the real component of the temperature for the zeroth homogeneous equations (section 3.5.1) - ^ X X J[cU,„,m,k][¥U),m0R,„ + < - a ' ) , m ^ „ ] d a ' • ft m=l n=l Truncating the number of vertical modes at N , and the wavenumber spectrum at A gives _ l y y frc1. JIV ., eR + V r R e1, Ida'. / i / . J |_ a,a ,n,m,k J|_ T (a-a ),m a ,n T (a-a),m"a,n J ft m=l n=l -A Consider IfC1 . Jwl 0 R +w R 61, Ida' J |_ a,a ,n,m,k J|_Y (a-a ),m a ,n Y (a-a ),m a ,n J -A Separating the integral over positive and negative domains, changing the variable of integration over the negative values, reversing the limits of integration, and using (3.28) this becomes ][cJ,J[sgn(a-a')V[a_, lmeRn + ¥ R_ a , l me:. n]da' +j[c;_, )][VU),m0'„ " K ^ e ^ d a ' 0 0 T E R M (1) T E R M (2) (3.50) where: C^.) = C a a , n m k , Cj_ a l ) = C a ( _ a . ) n m k , sgn(a - a') denotes the sign of ( a - a ' ) , and la - a' I denotes the absolute value of (a - a'). Chapter 3: Modified Rayleigh-Benard Convection ... 88 3.7.2 Numerical Solution Procedure Note that all of the analytical systems of governing equations may be written in the form | ^ = F(T,y) (3.51) where y = [ 0 R k , 0 a k , ¥ R , k > ¥ a , k ] T » a n d F contains both linear and nonlinear terms. It is the integral associated with the nonlinear terms which allows for the transfer of energy between various modes, and the majority of the computational effort is associated with determining the value of these integrals. Step 1: As illustrated in figure 3.4, the wavenumber spectrum is discretized into (J-1) equal intervals, involving J nodes, separated by a distance AJ = A / (J - 1 ) , where it is assumed that the coefficients associated with wavenumbers greater than A , are identically zero. Thus A represents the upper wavenumber limit. a=0 A j=l 2 3 4 J-2 J-1 J Figure 3.4: Discretization of the wavenumber spectrum. Since the number of vertical modes under consideration is truncated at N , the discretization of the governing equations, results in 4 * N * J equations to be numerically evaluated at each time step. Chapter 3: Modified Rayleigh-Benard Convection ... 89 Step 2: Initial values for the real and imaginary components of the streamfunction and temperature are prescribed. Typically all but the first vertical mode is assigned a null value. The three classes of initial conditions used and the values of the associated parameter sets are presented in detail in section 3.8.1. Step 3: Using the known values of the variables, the RHS of (3.51) may be evaluated explicitly at each time (or intermediate time) step. In order to solve the integrals associated with the function F of (3.51) a number of numerical methods were investigated including: Romberg integration with successive increases in the number of terms used to evaluate the integral based on a specified error tolerance, trapezoidal-type rule using all nodes with polynomial interpolation for the midpoint values, direct application of the trapezoidal rule using the nodal values, and finally the extended Simpson's rule (Press et al, 1992). The primary difficulty is that of computational speed as a result of the refinement process associated with the Romberg integration, and the number of calls to an interpolation subroutine when using a trapezoidal-type method combined with an mid-point interpolation scheme. In the case of the direct application of trapezoidal method, accumulated errors proved this technique impractical. As mentioned in the introduction of section 3.7, the most successful integration scheme, in terms of numerical accuracy and speed resulted from the use of the extended Simpson's rule, where an integral of the form A J[K(a')]da' o may be replaced by a summation over the nodal values Kj = K(a' i ) where a', = 0 and a', = A (figure 3.4) of the form Chapter 3: Modified Rayleigh-Benard Convection ... 90 / [ K ^ J d a ^ A J 1 4 2 4 2 4 1 - K , + - K 2 + - K , + - K . + . . . + - K . + - K , . + - K , .3 3 2 3 3 3 4 3 J 3 J~' 3 J + 0(J~ 4). (3.52) For example, consider an integral associated with the function F of (3.51), which is of the form T E R M (1) of (3.50), i.e. £ X j F m ( a - a ' ) G n ( a ' ) d a ' - (3-53) m=l n=l o The kernels of the integrals contained in both the homogenous and forced equations were evaluated by considering the sum of the products of the associated variables at each of the nodal values, for each of the vertical modes. Thus the kernel is defined as N N K = X X F m ( a - a ' ) G n ( a ' ) (3-54) m=l n=l where F m and G n represent the Fourier coefficients. Simplifying the notation by assigning the value of the node corresponding to (a-a ' ) with the index ' i j ' and the node associated with a' as ' j ' , (3.54) may be written N N K j =XXFiJ,mGj,n • m=l n=l Then expression (3.53) may be approximated using (3.51). Step 4: As with the evaluation of the integrals associated with the nonlinear terms, the forward time evolution of the solutions were evaluated using a number of methods. The adaptive, fifth-order embedded Runge-Kutta formulation using Cash-Karp parameters (Press et al., 1992), proved unsuccessful primarily as a result of the number of calls to the DERIVS subroutine which evaluates the RHS of (3.51) and involves computing the value of the nonlinear terms. Numerical Chapter 3: Modified Rayleigh-Benard Convection ... 91 tests conducted (with results presented in section 3.8) were done using both a finite-step version of the fifth-order Runge-Kutta method, as well as the original fourth-order Runge-Kutta scheme. Step 5: The final step is to ensure that the solution is converging towards steady-state. This is most readily done by considering the wavenumber development of the 'oldest' roll over the time period under consideration. In order to determine whether or not the solution has attained its steady-state value, the initial profile was integrated to some finite time (the length of which was determined in part by numerical limitations) with a series of up to 20 intermediate profiles output to file. The intermediate profiles were then Fourier transformed to mapped space, and the wavenumber of the 'oldest' roll is determined as a function of time along the half-way point between the two plates. Plotting the resulting time series, the relative stability of the wavenumber selection may be determined. The solution is mapped back to physical space and the wavenumber of the oldest roll is calculated in order to compare roll spacing with cloud street observations. The dominant scale associated with clear air roll formation is best determined in wavenumber space as data collected by aircraft is typically presented in spectral form. Results should not differ significantly, particularly if the preferred scale in wavenumber space is well defined. If not, a preferred range of wavenumbers may be estimated. Chapter 3: Modified Rayleigh-Benard Convection ... 92 3.8 Numerical Results Although the zeroth order equations will not provide insight into the question of the influence of the upper gravity wave field on scale selection as modeled, they do provide an opportunity to investigate the role of the nonlinear terms alone as a mechanism by which large aspect ratio rolls may be obtained. These results may be used to compare directly with laboratory experiments as the proposed mechanism for the formation of large aspect ratio rolls involving the coupling to an overlying gravity wave field, is obviously not plausible in this case. Comparison of results with those obtained by Getling (1983) provide a good test of the performance of the current model. Numerical results obtained from the first order equations (section 3.8.3) are corroborated by linear results presented in section 3.6. Tests based within the conductive regime demonstrate the instability of the conductive solution for all nonzero values of the Grashof number. Although computationally very 'fast' when compared to the numerically intensive truncated equations (section 3.8.4), the shortcomings of the perturbation analysis (to first order) are demonstrated via comparison with the results from the truncated analysis (section 3.8.5). 3.8.1 Spectral Initialization In general, nonlinear systems are highly dependent on the initial conditions. Therefore, even though the governing equations may accurately describe the time evolution of the flow, the dependence of the nonlinear system on exact specification of the initial conditions represents the main obstacle to long term predictability. Since in practice the initial state may be defined only to within some degree of uncertainty, sensitivity tests play an important role in developing confidence in the results obtained from a numerical model of a nonlinear system. Chapter 3: Modified Rayleigh-Benard Convection ... 93 Numerical simulations require not only specification of the initial form of the temperature and streamfunction perturbation, but also the amplitude (e) and wavenumber (b) associated with the warped upper surface. We shall see in section 3.8 that the conclusions relating to wavenumber selection are insensitive to the choice of b, and a range of wavenumbers (all of which are less than ac) are considered. Al l numerical simulations are conducted using a value of 8 of 0.01 unless otherwise specified. This value is arbitrary and assumes that the amplitude of the undulations are 1% of the depth of the boundary layer and is believed to be a conservative estimate of the possible effect of the gravity waves on the inversion. We shall see (section 3.8.4) that for a broad range of values of the amplitude of the warping, the results associated with wavenumber selection are robust as the conclusions drawn from the analysis are shown to be insensitive to the choice of 8. However, the magnitude of e will influence the energetics of the roll circulations and therefore net heat transport through the domain (sections 3.8.4 and 3.9). The behavior of the solution for the streamfunction and temperature is investigated using three types of initial spectra in wavenumber space. Class I initial conditions (Getling, 1983) assumes that energy is available to a broad spectrum of horizontal wavenumbers (as opposed to randomly distributed energy) covering the unstable range of wavenumbers as predicted by linear theory (Rayleigh, 1918). Assuming that energy exists over a broad spectrum of wavenumbers allows the flow 'equal access' to a wide range of horizontal length scales. Class I initial conditions represent an initially localized roll and temperature perturbation in physical space. Class II initial conditions (Getling, 1983) assumes the initial state is one of pre-existing, horizontally periodic temperature and flow perturbation, and is equivalent to the distribution of energy in a narrow band in wavenumber space. Both of these class types attempt to simulate probable initial scenarios Chapter 3: Modified Rayleigh-Benard Convection ... 94 which may be reproduced in laboratory experiments. Of the two classes proposed by Getling, it is reasonable to assume that the class II initial conditions would more accurately represent the atmospheric scenario associated with the initial generation of the gravity wave field. Finally class III initial conditions are used to study the temporal development of the first order and truncated equations, and simply represent the warped temperature field with no superimposed roll or temperature perturbation. Class I The real and imaginary components of the temperature and streamfunction are assigned initial values corresponding to ¥ a R , = - K v (0<a<a 2 ) - K v ( a 3 - a ) / ( a 3 - a 2 ) ( a 2 <a<a 3 ) 0 (a3 < a < A) f K v a/a, K v K v ( a 3 -a) / (a 3 - a 2 ) 0 (0<a<a,) (a! <a<a 2 ) (a 2 <a<a 3 ) (a 3 < a $ A) with \ j / R k = \|/ a k = 0 for k greater than one. The real and imaginary components of the temperature take the same form, with K v replaced by (-KT) for the real component, and K v replaced by K T for the imaginary component. The energy is assumed to be contained entirely by the first vertical mode, as all coefficients associated with the higher vertical mode numbers are assigned a null value. See figure 3.5(a) and 3.5(b) for the initial conditions in (%,TI) space. Chapter 3: Modified Rayleigh-Benard Convection 95 -10 -5 0 5 10 Figure 3.5(a): Class I initial conditions in mapped space: Streamfunction. Contour intervals: -0.08, -0.06, -0.04, -0.02, 0.0, 0.02, 0.04. Vertical lines denote the zero contour. Alternating rolls have positive and negative values of the streamfunction with the dot • denoting negative values of the streamfunction. Parameter set 'a'. o -10 -5 0 5 10 Figure 3.5(b): Class I initial conditions in mapped space: Temperature. Contour intervals: -1, -0.75, -0.5, -0.25, 0.0, 0.25, 0.5, 0.75, 1. Note that the contours increase in value from top to bottom. Parameter set 'a'. Chapter 3: Modified Rayleigh-Benard Convection ... 96 The numerical tests which were conducted using the class I initial profile had values of the prescribed parameters as shown in table 3.1. For abbreviation the notation 'Cla ' (for example) denotes that a particular numerical test was run using class I initial conditions with set 'a' parameter values. Parameter set 'a' and 'g' were chosen to span the linear stability range, sets 'b' and 'c' are used to demonstrate the effects of the choice of amplitude of the imposed perturbations on wavenumber selection, while parameter sets 'd' and 'e' are used to investigate the sensitivity of the obtained solutions to the width of the initial, localized roll. Finally, set ' f is used in the analysis of the first order and truncated equations, in order to compare the amplitude of forced and unforced wavenumber development. Set K V KT ai a 2 a3 a 0.050 0.050 1.25 5.0 5.4 b 0.050 0.005 1.25 5.0 5.4 c 0.005 0.050 1.25 5.0 5.4 d 0.050 0.050 1.25 3.0 5.4 e 0.050 0.050 1.25 4.0 6.4 f 0.050 0.050 1.25 6.0 6.4 g 0.050 0.050 1.25 9.0 9.4 Table 3.1: Class I parameter values. Class R Following Getling (1983), the real and imaginary components of the streamfunction and temperature are assigned the initial values 4sV7t ( a -a ) Chapter 3: Modified Rayleigh-Benard Convection ... 97 with k = \|/ a k = 9 a k = 6 a k =0 for k greater than one. The variables a , S 0 and s are external constants. In the limit as s —> 0, the class II initial profiles represent a delta function in wavenumber space. See figures 3.6(a) and 3.6(b) for the initial profiles of the streamfunction and temperature in (£, r\) space. Table 3.2 lists the parameter values used in the numerical test. As with the class I parameter sets, sets 'b' and 'c' where chosen in order to investigate the effect of the external parameters on scale selection, while'd' and 'e' are used to investigate the effect of the initial roll width. Parameter set ' f is used in the initialization of the truncated equations with values of the parameters chosen such that the amplitude of the forcing and the maximum amplitude associated with wavenumbers other than the imposed wavenumber and it's harmonics are comparable for the time under consideration. (See section 3.8.4.) Set a S s b d e c a f 2.2 0.355 0.400 2.2 0.355 0.100 2.2 0.100 0.400 1.0 0.355 0.400 3.0 0.355 0.400 2.2 0.100 0.100 2.4 0.005 0.100 Table 3.4: Class II parameter values. Chapter 3: Modified Rayleigh-Benard Convection 98 Figure 3.6(a): Class II initial conditions in mapped space: Streamfunction. Contour intervals: -0.09, -0.06, -0.03, 0.0, 0.03, 0.06, 0.09. Vertical lines denote the zero contour. The dot • denotes negative values of the streamfunction, with adjacent rolls alternating in sign. Parameter set 'a'. Figure 3.6(b): Class II initial conditions in mapped space: Temperature. Contour intervals: -1, -0.75, -0.5, -0.25, 0.0, 0.25, 0.5, 0.75, 1. Note that the contours increase in value from top to bottom. Parameter set 'a'. Chapter 3: Modified Rayleigh-Benard Convection ... 99 Class in In addition to the two class types suggested by Getling, a third class type was used to initialize the first order equations and the truncated equations. For the case of the warped upper surface, the class III initialization consists simply of the warped temperature field with no initial streamfunction or temperature perturbation. This initial condition is highly idealized since there is no initial energy made available to scales other than those associated with 'b' and its harmonics. Class III initial conditions are depicted in figures 3.7(a) and 3.7(b) in mapped space and physical space respectively. 3.8.2 Results from the Zeroth Order Equations As discussed in section 3.3.3, the zeroth order equations govern the nonlinear dynamics associated with the classical R-B convection problem and has been studied by a number of researchers (see section 2.2). In particular, the current study follows the analytical development and numerical solution techniques of Getling (1983) and comparison of model results with those of Getling provided a good test case for verification of model performance. The objectives of the numerical tests were to • determine the behavior of the flow for supercritical values of the Grashof number over a range of Prandtl number (0.1, 1., 10), in order to determine the trend of the preferred scale; (Note that the values of G and a were chosen in order compare with results of Getling (1983).) • conduct a sensitivity study in order to evaluate the influence of external model parameters on results. • investigate the sensitivity of results obtained to initial conditions; Chapter 3: Modified Rayleigh-Benard Convection 100 -10 J 0 0 10 Figure 3.7(a): Class III initial conditions in mapped space: Temperature. Contour intervals: -1, -0.75, -0.5, -0.25, 0.0, 0.25, 0.5, 0.75, 1. Note that the contours increase in value from top to bottom. Parameter values: b = 1, e = 0.025 . Note that the amplitude of the forcing was chosen merely for illustration purposes. -10 0 X -» (l + ecos(bX)) T Z 0 10 Figure 3.7(b): Class III conditions in physical space: Temperature. Contour intervals: -1, -0.75, -0.5, -0.25, 0.0, 0.25, 0.5, 0.75, 1. Note that the contours increase in value from top to bottom. Parameter values as in figure (3.7a). Chapter 3: Modified Rayleigh-Benard Convection ... 101 Summary of Model Results • Sensitivity Study: 0 Vertical mode (N) and horizontal wavenumber (A) cutoff: Results suggest that for all values of the parameters tested (with the exception of a Prandtl number of 0.1), for the length of the tests conducted, the choice of vertical modes and horizontal wavenumber cutoff were adequate to ensure confidence in the wavenumber development. However, results of the small Prandtl number (0.1) tests which were conducted for 1.52 times the critical value of the Grashof number were found to be very dependent on the value of the vertical mode and horizontal wavenumber truncation parameters. Truncation effects have been investigated by Fujimura and Mizushima (1988) in association with a Fourier series expansion for the problem of convection between vertical parallel plates. Their results indicated that qualitatively different solutions were obtained which depended on the number of terms used in the Fourier series expansion. ° Class initialization parameter values: Class I: Results suggest that as long as the initial breadth of the spectrum includes the (ultimately) chosen scale, flow development adjusts quickly. An order of magnitude variation of the initial amplitude of the streamfunction (or temperature) perturbation showed only slight sensitivity (in particular to changes in the thermal amplitude) to the choice of the parameter K v , and K T . Class II: Although the results are insensitive to the amplitude of the initial perturbation, class I results are best summarized by suggesting that care be taken when selecting the bandwidth and initial central wavenumber if the approximate value of the final preferred wavenumber is not known a priori. ° Nodal resolution: The stability of the solutions were fairly insensitive to the number of nodes (J) chosen for the duration of the tests conducted. Care must be taken to ensure that errors Chapter 3: Modified Rayleigh-Benard Convection ... 102 associated with the integration over wavenumber space are controlled. Getling (1983) discusses the development of 'peaks' which result from errors associated with the use of Simpson's rule when evaluating the nonlinear integrals. As a result of the limitation of Simpson's method when applied to rapidly oscillating functions, errors develop as the oscillations of the spectra approaches nodal resolution. Contamination of spectra as a result of numerical errors will be further discussed in section 3.9 which outlines numerical modeling efforts at values of the Grashof number associated with atmospherically observed phenomena. • Initial Conditions: For all cases investigated, the roll development (in mapped space) was found to be highly dependent on the form of the initial state specified. As the flow field develops, the existence of neighboring rolls inhibits rapid cell expansion or contraction. In association with the form of the initial class II conditions, results indicate a very slow transition to the preferred scale if the bandwidth is too narrow. Since little or no energy may be contained at the preferred scale, nonlinear transfer of energy and subsequent wavenumber development will be inhibited. • In agreement with previous studies, results suggest that a two-dimensional mechanism is not able to account for roll development with scales larger than the critical wavenumber for Prandtl numbers of 1 (Getling, 1983; Rothermal and Agee, 1986) and 10 (Getling, 1983) for G / G c = 1.004 and G / G c = 1.52. However, there is a noted shift to smaller wavenumbers (i.e. larger scales) for a Prandtl number of 0.1 as found by Getling (1983) and Chang and Shirer (1984). Chapter 3: Modified Rayleigh-Benard Convection ... 103 3.8.3 Results from the First Order Equations Beyond the analytical results presented in section 3.6.2, numerical tests based on the first order equations provide a more detailed look at the effects of the corrugated upper surface on flow dynamics. However, the noted linearity of the first order equations (section 3.3.4) has lead to limited numerical testing. A comparison of results with those obtained from the solution to the truncated equations is presented as part of section 3.8.5. Recalling the first order equations given by equations (3.35) and (3.36), as well as the assumptions associated with the form of the zeroth order solutions outlined in section 3.6.2, results for subcritical ( G < G C ) and supercritical ( G > G C ) values of the Grashof number are presented. The values of the numerical test parameters are given in table 3.3. Test a G T f i N A nodes step size method F120 1 1 1 1 6 121 lxlO" 2 rk4. Cla F l 1 1 2 1 6 121 lxlO" 2 rk4, Cla F100 1 1 2 3 6 121 l x l O 2 rk4, Cla F2 1 250 1 1 6 121 lxlO" 2 rk4, Cla F200 1 250 2 1 6 121 lxlO" 2 rk4, Cla F210 1 250 4 1 6 121 lxlO- 2 rk4, Cla F110 1 250 2 3 6 121 lxlO" 2 rk4, Cla F331 1 330 2 1 6 121 lx lO ' 2 rk4, CIII F332 1 330 4 1 6 121 l x l O 2 rk4, CIII F333 1 330 6 1 6 121 l x l O 2 rk4, CIII F330 1 330 12 1 6 121 lxlO" 2 rk4, CIII F340 1 330 12 3 6 121 lxlO" 2 rk4, CIE Table 3.3: First order equations. Numerical test parameters. Chapter 3: Modified Rayleigh-Benard Convection ... 104 Conductive Base State The case G < G c is of particular interest since the numerical model verifies analytical results which indicate that a secondary circulation is induced for all (nonzero) values of the Grashof number. Figure 3.8 depicts the solution of the first order equations in wavenumber space for test F l using Cla initial conditions. The value of the parameters correspond to G = l , a = 1, b = 1 and N = l , with A R = A 1 =0 , B R = B r =0 and a final nondimensional time of 2.0 (or approximately 80 minutes). Although the amplitude of the forcing does not appear explicitly in the first order equations, for the purpose of mapping the solution back to physical space, a value of 0.01 was used. For the vector of unknowns given by [0R,9J,\|/R,\|/|]T, analytical results (section 3.6.2) predict a forcing of f, = [-2,0,0,2]T and a solution vector with component values of [-0.3748, 0.0, 0.0,-0.0370]T, which was corroborated with numerical predicted values. (See appendix E for details of comparison between the linear analysis and numerical results). The streamfunction and temperature are illustrated in mapped and physical space in figures 3.9(a) and 3.9(b) respectively. A comparison of analytical and numerical test results for a Grashof number of 250 (with o" = 1) was also conducted again showing excellent agreement with details given in appendix E. Chapter 3: Modified Rayleigh-Benard Convection 105 ISTFI o T3 H S 0.04 0.035 h 0.03 h 0.025 0.02 0.015 h 0.01 0.005 h 0 ITempI J_ _L 0 1 2 3 4 wavenumber •o 3 E C3 1 2 3 . 4 5 wavenumber Figure 3.8: Results from the first order equations in wavenumber space. Conductive regime. Magnitude of the streamfunction ISTFI and temperature ITempI. Parameter values: G = 1, o" = 1, b = 1 and N = 1. Chapter 3: Modified Rayleigh-Benard Convection 106 -5 0 5 10 X -» Figure 3.9(a): Solution to the first order equations. Streamfunction. Top: Streamfunction in mapped space. Bottom: Streamfunction in physical space. Contour intervals: -2.5c"4, -1.25e"4, 0, 1.25e"4, 2.5e"4. Vertical lines denote the zero contour and the dot • denotes negative values of the streamfunction. Parameter values: G = 1, <J = 1, b = 1 and N = 1. Chapter 3: Modified Rayleigh-Benard Convection 107 T -10 -5 _1 0 0 10 (l + ecos(bX)) t Z 0 -10 0 X -> 10 Figure 3.9(b): Solution to the first order equations. Temperature. Top: Temperature in mapped space. Bottom: Temperature in physical space. Contour intervals: -1, -0.75, -0.5, -0.25, 0.0, 0.25, 0.5, 0.75, 1. Note that the contours increase in value from top to bottom. Note that the perturbation is so small that it does not appear in temperature field. Parameter values as in figure 3.9(a). Chapter 3: Modified Rayleigh-Benard Convection ... 108 Convective Base State Figure 3.10 depicts the solution of the first order equations in wavenumber space for a G / G c = 1.004 using CIII initial conditions with a = 1, b = 1, N = 1. The amplitudes of the components of the zeroth order variables used to evaluate the forcing functions h ( i | / ° k , 6 ° k ) and q ( ¥ a k ' ° a k ) ' w e r e t a ken from the results of the zeroth order equations and are assigned the values A R = A 1 =0.04525 and B R = ( -B 1 ) = 0.1511. For the vector of unknowns given by [0R,0j,\|/R,V|/J]T , analytical results of section (3.6.2) predict forcings of "-4 T " 0.4160" T ' 0.4658" 0 0.4160 0.4658 f, = f2 = u = 1 0 z 28.305 3 16.670 1320_ -28.305 -16.670 and corresponding solution vectors "-4.986 " e b 0 0 -25.100_ A R D ( a „ - b ) " 0.4014" ~A R U ( a „ + b ) " 0.5086 A 1 ° ( a „ - b ) 0.4014 A1 ° ( a „ + b ) 0.5086 < 0 - b ) -1.5035 V ( a 0 + b ) -1.4682 <.-•»_ 1.5035 1.4682 . V ( . 0 + b ) . which agree well with numerical predicted values. (See appendix E for details). The solution vectors obtained over 97% of their final value by timet = 2 , 99.9% byx = 4 , and approximately 100% (to three decimal places) by T = 6. Plotted in the top two figures of 3.10, are the magnitude of the streamfunction ISTFI and the magnitude of the temperature ITempI of the zeroth order variables used to evaluate the forcing functions. The lower two figures are the resulting solutions of the first order equations. Note that Chapter 3: Modified Rayleigh-Benard Convection 109 2 3 4 5 wavenumber -a "H. £ CO 0.07 0.06 h 0.05 0.04 h 0.03 0.02 -0.01 -0 _L ZOE: Itempl i 1 r _L 0 1 2 3 4 5 wavenumber 2 3 4 5 wavenumber ITempI •o 3 a* B CO 1 2 3 4 - 5 wavenumber Figure 3.10: Results from the first order equations in wavenumber space. Convective regime. Top: Zeroth order solution used to force the first order equations. Shown: Magnitude of the streamfunction and temperature. Bottom left to right: Solution of the first order equations; magnitude of the streamfunction, magnitude of the temperature. G = 330 , a = 1, b = 1 and N = l . Chapter 3: Modified Rayleigh-Benard Convection ... 110 the dominant scale is associated with the imposed wavenumber (b = 1) and results from the delta function; i.e. term (1), equation (3.35). The two smaller peaks at a = 1.2 and a = 3.2, result from the shifting left and right in wavenumber space of the zeroth order solution, which was centered about a 0 = 2.2. Recall however, that the complete solution as predicted by the perturbation analysis is the sum of the zeroth order solution and e times the first order solution. Results from these two sets of equations predict that as a minimum, multiple scaling will result from the warping of the inversion. First order equation results suggest that which mode will form the dominant one is dependent on whether the imposed scale or the preferred scale associated with flow development as governed by the zeroth order equations contains more kinetic energy. These results will be shown to be misleading when compared to the results of the truncated equations presented in the next section. 3.8.4 Results from the Truncated Equations Recall that the truncated equations govern the nonlinear flow dynamics as influenced by the thermal forcing associated with the rigid, isothermal, perfectly conducting, corrugated upper surface, and that the equations were obtained via expansion of the coordinate mapping to first order in e. Of particular interest is the form of the solution as a function of the imposed wavenumber 'b'. Due to the increased computational requirements (c.f. with evaluating the first order equations) tests were conducted for o = 1 only, with the numerical test parameters given in table 3.4. Chapter 3: Modified Rayleigh-Benard Convection ... Test a G b T f N A nodes step size method C l I 1 1.0 4. 1 6. 121 5xl0" 3 finite. Cla C100 I 1 1.0 4. 3 6. 121 5 x l 0 3 finite, Cla C2 I 250 1.0 8. 1 6. 121 5xl0" 3 finite, Cla C200 I 250 1.0 8. 3 6. 121 5xl0" 3 finite, Cla C30 I 330 1.00 12. 1 6. 121 5xl0" 3 rk4, CIII C31 I 330 1.50 12. 1 6. 121 5xlO"3 rk4, CIII C32 I 330 1.75 12. 1 6. 121 5 x l 0 3 rk4, CIII C33 ] I 330 2.00 12. 1 6. 121 5xl0" 3 rk4, CIII C34 1 [ 330 1.00 12. 3 6. 121 5xl0" 3 rk4, CIII C35 1 I 330 1.50 .12. 3 6. 121 5xl0" 3 rk4, CIII C36 ] I 330 1.75 12. 3 6. 121 5xlO"3 rk4, CIII C37 ] I 330 2.00 12. 3 6. 121 5xl0"3 rk4, CIII C400 1 I 330 1.00 12. 1 6. 121 lxlO" 2 rk4, C l l f C410 1 [ 330 1.50 12. 1 6. 121 lxlO" 2 rk4, C l l f C420 ] I 330 1.75 12. 1 6. 121 l x l O 2 rk4, C l l f C430 1 [ 330 2.00 12. 1 6. 121 lxlO" 2 rk4, C l l f C440 ] [ 330 1.00 12. 3 6. 121 lxlO" 2 rk4, C l l f C450 ] I 330 1.50 12. 3 6. 121 l x l O 2 rk4, C l l f C460 ] I 330 1.75 12. 3 6. 121 l x l O 2 rk4, C l l f C470 ] [ 330 2.00 12. 3 6. 121 lxlO" 2 rk4, C l l f C50 1 I 500 1.00 4. 1 10. 201 lxlO" 2 finite, CHI C51 1 I 500 1.50 4. 1 10. 201 l x l O 2 finite, C m C52 1 [ 500 1.75 4. .1 10. 201 lxlO" 2 finite, C m C53 [ 500 2.00 4. 1 10. 201 lxlO" 2 finite, C m C54 1 I 500 1.00 4. 3 10. 201 lxlO" 2 finite, CIII C55 I 500 1.50 4. 3 10. 201 lxlO" 2 finite, CIII C56 1 [ 500 1.75 4. 3 10. 201 l x l O 2 finite, CHI C57 1 [ 500 2.00 4. 3 10. 201 lxlO" 2 finite, CIH C60 ] I 500 1.00 4 1 10 201 lxlO" 2 finite, CHg C61 I 500 1.50 4 1 10 201 l x l O 2 finite, Cl lg C62 I 500 1.75 4 1 10 201 lxlO" 2 finite, CHg C63 I 500 2.00 4 1 10 201 lxlO" 2 finite, CHg C64 : 5 500 1.00 4 3 10 201 l x l O 2 finite, CHg Table 3.4: Truncated equations. Numerical test parameters. Chapter 3: Modified Rayleigh-Benard Convection ... 112 Conductive Regime For G < G c , numerical results agree with those obtained from the linear analysis which indicated multiple scaling (section 3.6.3 and appendix E) associated with the imposed wavenumber 'b' and its harmonics. Figure 3.11 illustrates the final wavenumber spectra for the ISTFI for test C100 with G = 1 and b = 1, and for C200 with G = 250 and b = 1. Although the ISTFI associated with the forced circulation is small, it is nonzero at 'b' and the harmonics of b. Both tests were initialized using class I initial conditions which represents a localized roll in physical space. This initial disturbance quickly dissipated as the only scales of motion which are supported by the flow for G < G c are those associated with 'b' and its harmonics. Convective regime For G > G c , two classes of initial conditions were used: class III and class II. The objective of numerical simulations conducted using the class III initial conditions was to focus on flow development solely in response to a warping of the isotherms. The objective of test based on class II initial condition was to consider the effects of the warping on a pre-existing set of roll vortices. a) Class III initial conditions Recall that class III initial conditions (section 3.8.1) represent a warped temperature field with no superimposed perturbations. This profile is highly idealized since no energy is supplied to any wavenumber other than those associated with 'b' and its harmonics. Thus the null value associated with other scales will be very slow to develop and will only attain nontrivial values in association with energy transfer resulting from such terms as T E R M (3) of equation (3.35). Chapter 3: Modified Rayleigh-Benard Convection 113 ISTFI -o 3 Cu E 0.0004 0.0003 0.0002 o.oooi H J L 0 0 2 4 6 wavenumber -a ISTFI 0 2 4 6 wavenumber -o 4—1 p 5 le-07 8e-08 6e-08 4e-08 2e-08 0 0 2 4 6 wavenumber ISTFI •o 3 "S. e 0.0004 0.00035 0.0003 0.00025 0.0002 I-0.00015 h 0.0001 5e-05 h 0 0 1 _L 2 3 4 wavenumber Figure 3.11: Wavenumber spectra. Conductive regime. t Top left to right: ISTFI, and close-ups showing nonzero values of the streamfunction at the harmonics of the imposed wavenumber b = 1 for G = 1. Bottom: ISTFI for G = 250 and b = 1 .Other parameter values: o~=l, b = l , e = 0.01, T f = 4 with Cla initial conditions. Chapter 3: Modified Rayleigh-Benard Convection ... 114 In practice, class III initial conditions are not realizable since imperfections (such as surface nonuniformities) will make energy available to a wide range of scales. Nonetheless, it does allow for a focused look at the effects of the forcing. In particular, class III initial conditions introduce a horizontal length scale based on the imposed wavenumber 'b' and flow development is dominated by 'b' and its harmonics. Since the horizontal dependence of the upper surface is represented using a cosine function, flow development using class III initial conditions tends to isolate \ | / a k which is associated with sin(b^), and 0 a k which is associated with cos(b^); i.e. \ | / R k = 0 a k =0 (see figure 3.12). This is not unexpected since the temperature is then in phase with the upper warping and the streamfunction is 90 degrees out of phase resulting in the upward transport of heat from the warm lower surface to the cold upper surface. The value of the ISTFI for the first vertical mode resulting from linear and nonlinear numerical simulations with G / G c = 1.004 (i.e. G = 330), are presented in table 3.5. Similarly for G / G c = 1.52 (i.e. G = 500), linear and nonlinear results are outlined in table 3.6. Tests conducted with G / G c = 1.004 indicate only slight differences between linear and nonlinear results. This is not surprising since the value of the Grashof number is only marginally above the critical value. This suggests that flow dynamics may be adequately described by a linear model as the amplitude of the resulting circulation is small and thus nonlinear influences will be negligible. Chapter 3: Modified Rayleigh-Benard Convection 115 CD T3 3 Cu 0 2 4 wavenumber Im(STF) 0.1 0 -0.1 T3 13 -0.2 "CU E -0.3 -0.4 -0.5 2 4 wavenumber 0.05 -0.05 h 0 Re(Temp) T 2 4 wavenumber -a 3 •4—» "cu Im(Temp) T 2 4 wavenumber Figure 3.12: Wavenumber spectra obtained from the truncated equations. Convective regime. Top left to right: Real and imaginary components of the streamfunction. Bottom left to right: Real and imaginary components of the temperature. Parameter values: G = 330, N = 3, A = 6, b = 1, e = 0.01, T f = 12 and CIII initial conditions. Chapter 3: Modified Rayleigh-Benard Convection . 116 Linear Non Linear b < l m a x ISTFI ISTFI ISTFI Test ISTFI ISTFI ISTFI Test a=b a=2b a=3b a=b a=2b a=3b 1.00 2.00 2.78e_1 4.02c"1 5.82e'2 C30 2.786"1 4.05c 1 5.86e'2 C34 1.50 1.50 1.00 1.45c'1 1.43e"4 C31 1.00 1.45e"] 2.44e"3 C35 1.75 1.75 2.70 1.41e"' 1.27e"3 C32 2.70 1.426"1 1.31e"3 C36 2.00 2.00 10.44 2.72e"! - C33 10.36 2.7061 - C37 Table 3.5: Linear (N = 1) and nonlinear (N = 3) streamfunction amplitudes, first vertical mode. Parameter values: G = 330, e = 0.01, T f = 12 , and CIII initial conditions. Linear Non Linear b ^ m a x ISTFI ISTFI ISTFI Test ISTFI ISTFI ISTFI Test a=b a=2b a=3b a=b a=2b a=3b 1.00 2.00 2.01e2 4.38e3 1.86e3 C50 1.37c1 3.63e2 7.03 C54 1.50 3.00 1.7 le 3 8.66e3 1.85e2 C51 5.93e: 2.27e2 6.25 C55 1.75 1.75 1.37e4 1.97e3 2.16c1 C52 2.58e2 3.01 1.62 C56 2.00 2.00 1.53e5 5.58e3 - C53 3.65e2 3.36 - C57 Table 3.6: Linear (N = 1) and nonlinear (N = 3) streamfunction amplitude, first vertical mode. Parameter values: G = 500, e = 0.01, T f = 4 , and CIII. Note that for b = 1 linear results (C30) (slightly) under estimates the ISTFI associated with the first harmonic of 'b', i.e. at a = 2. Based on two-dimensional numerical simulations the preferred wavenumber (a p ) for classical R-B convection (zeroth order equations) is estimated at Chapter 3: Modified Rayleigh-Benard Convection ... 117 2.2 for G / G c = 1.004 (and 2.4 for G / G c = 1.52). Results from table 3.6 illustrate that the influence of the warping on flow dynamics increases as the imposed wavenumber approaches a p . However it will the harmonic closest in wavenumber space to a p which dominates. Also note that although both C33 and C37 indicate the dominance of the wavenumber a = 2, there is a difference of two orders of magnitude in the value of the ISTFI depending on whether the dominant wavenumber is associated with the imposed wavenumber (C34) or the first harmonic (C37). Comparison of linear and nonlinear results for G / G c = 1.52 indicate that the numerical simulations based on a linearized version of the truncated equations, results in an over estimation of the magnitude of the flow variables. This is not surprising since the inherent nonlinearity of the system will transfer energy from initially dominant scales (as determined by linear theory) to other scales of motion. Recall that the objective of the study was to investigate gravity wave-boundary layer interactions as a possible mechanism for the formation of large aspect ratio (i.e. > 4) rolls. Results given in table 3.5 ( G / G c = 1.004) and table 3.6 ( G / G c = 1.52) indicate that although scales of motion larger than the wavelength observed at the critical Grashof number (X c ~ 2.8) may be obtained in response to the warping of the upper surface, the dominant scale will be associated with the harmonic closest (in wavenumber space) to the preferred scale of classical R-B convection, and thus may not be associated with the forced scale. This suggests that the dominant scale 'a' will lie in the range ( 2 / 3 < a / a p <4/3) where ap is the preferred wavenumber as Chapter 3: Modified Rayleigh-Benard Convection ... 118 selected by convection between flat plates (i.e. classical R-B convection). Based on this range of values a maximum increase in roll scales (i.e. X) of 50% is predicted. In order to test the robustness of these results, numerical simulations were conducted with G / G c = 1.52 for values of the amplitude of the warping given by 0.001, 0.01, and 0.1. Other parameter values include an imposed wavenumber of 1 and a final nondimensioanlized time of 4. Results depicted in figure 3.13 indicate the dominance of the first harmonic associated with the wavenumber 2, in all three cases. Results from the nonlinear test C54 (b = 1) and C55 (b = 1.5) are presented in figure 3.14. The top two figures depict the ISTFI and dominance of the second harmonic is evident in both cases. The middle and lower panel presents the streamfunction in physical space for both cases in order to compare observed roll geometry. Note that the number of rolls contained in one imposed wavelength is four for the middle and lower panels as expected (i.e. (X = 7t / b) where X is the observed wavelength). Application of normal mode analysis under a linear framework suggests that the fastest growing mode will dominate. However, figure 3.15 illustrates that the (initially) fastest growing mode (associated with b = 1) only dominates during the early stages of flow development. In particular, for both C34 ( G / G c = 1.004, b = 1) and C54 ( G / G c = 1.52, b = 1) the preferred scale is not associated with 'b', but rather the first harmonic of 'b' i.e. at a = 2b = 2 (tables 3.5 and 3.6) and adjustment of the flow from a wavenumber of 1 to 2, is evident in both figures. Noting the difference in the time scales, the preferred harmonic manifests itself more rapidly for the higher Grashof number case. Chapter 3: Modified Rayleigh-Benard Convection 119 3 6 400 300 h 200 h 0 1 2 3 4 5 wavenumber ISTFI ISTFI 1 2 3 4 5 wavenumber 400 300 h 3 ~ 200 h 0 1 2 3 4 5 wavenumber Figure 3.13: Effect of e on the solution obtained from the truncated equations. Top: Magnitude of the streamfunction ISTFI with an £ value of 0.01. Bottom left to right: ISTFI with e values of 0.001 and 0.1. Parameter values: G = 500 , N = 3, A = 6, T f = 4 , and Cffl. Chapter 3: Modified Rayleigh-Benard Convection 120 400 i 1 1 1 350 - -300 - -o 250 -3 200 — B 150 — -100 -50 0 . 1 I i 4 6 wavenumber 250 i 1 1 1 200 - -o ~o 150 3 "5. E a 100 -50 -0 1 10 4 6 wavenumber 10 (nl (nl (nl 1^1 (ni (nl (nl (nl (old U U i y y y U y y U u y ! -10 -5 (l + ecos(bX)) T Z 0 0 X - » A A A A A A A J w A A 10 _0 + £cos(bX)) t Z 0 -10 0 X - » 10 Figure 3.14: Comparison of results from truncated equations based on imposed wavenumber b, in physical space. Convective regime. Top left to right: ISTFI with b = 1 and b = 1.5. Middle: STF in physical space for b = 1 with maximum amplitude at a = 2. Contoured at intervals -2.5, -1.25, 0.0, 1.25, and 2.5. Bottom: STF in physical space for b = 1.5 with maximum amplitude at a = 3. Contoured at intervals -1.8, -0.9, 0.0, 0.9, and 1.8. The dot • denotes negatively contoured values of the streamfunction. Note the multiple-scaling. Parameter values: G = 500, N = 3, A = 10, e = 0.01, T f = 4 , and CIII. Chapter 3: Modified Rayleigh-Benard Convection 121 CU JO a 3 C CD > nondimensionalized time 2 3 nondimensionalized time Figure 3.15: 'Oldest' roll wavenumber development Top: Test C34 with G = 330, N = 3, b = l , 8 = 0.01, and T f = 12 . Bottom: Test C54 with G = 500, N = 3, b = 1, 8 = 0.01, and T f = 4 . Note the change in preferred scale. Chapter 3: Modified Rayleigh-Benard Convection ... 122 b) Class II initial conditions As previously noted, class III initial conditions are highly idealized. Therefore, in order to investigate the robustness of the results presented; namely the anticipated range of the dominant wavenumber, a second initialization class was used; CII with parameter sets T and 'g ' . Recall that class II initial conditions are characterized by a pre-existing set of rolls of a specified wavenumber in physical space (see section 3.8.1). Numerical tests were designed to identify the effect of the warped upper lid on a 'uniform' roll field. With respect to the selection of the dominant harmonic, results presented in table 3.7 and 3.8 for G / G c = 1.004, and in tables 3.9 and 3.10 for G / G c = 1.52, agree with those presented using CIII initial conditions (cf. tables 3.5 and 3.6); namely the harmonic closest to the preferred scale of the flat-plate case, dominates. b amax &o ISTFI at ISTFI at ISTFI at ISTFI at a=b a=2b a=3b a=a0 Test 1.00 2.00 2.2 1.50 1.50 2.2 1.75 1.75 2.2 2.00 2.00 2.2 2.78e"' 4.02c 1 5.82e"2 3.78c"1 1.00 1.45e"' 2.43e"3 2.98c"1 2.70 1.41e"' 8.02e"3 2.83c"1 10.44 2.72c"1 - 2.73c"1 C400 C410 C420 C430 Table 3.7: Linear (N = 1) streamfunction amplitudes, first vertical mode. Parameter values: G = 330 , 8 = 0.01, T f = 12 , and CITf. Chapter 3: Modified Rayleigh-Benard Convection . 123 b a m a x ao ISTFI at ISTFI at ISTFI at ISTFI at a=b a=2b a=3b a=a0 Test 1.00 2.00 2.2 1.50 1.50 2.2 1.75 1.75 2.2 2.00 2.00 2.2 2.78c"1 4.05c"1 5.87e"2 3.78c"1 1.00 1.456"1 2.44e"3 2.98c"1 2.70 1.42c"1 1.31e"3 2.81c"1 10.36 2.70e"' - 2.64c"1 C440 C450 C460 C470 Table 3.8: Nonlinear (N = 3) streamfunction amplitudes, first vertical mode. Parameter values: G = 330 , e = 0.01, T f = 12 , and Cflf. b ^max a<, ISTFI at ISTFI at ISTFI at ISTFI at Test a=b a=2b a=3b a=a0 1.00 2.4 2.4 2.01e2 4.38e2 1.86e3 1.85e4 C60 1.50 2.4 2.4 1.71e3 8.66e3 1.85e2 1.73e4 C61 1.75 2.4 2.4 1.37e4 1.97e3 2.16c1 1.70e4 C62 2.00 2.0 2.4 1.53e5 5.58e3 - 1.69e4 C63 Table 3.9: Linear (N = 1) streamfunction amplitudes, first vertical mode. Parameter values: G = 500, e = 0.01, T f = 4 , and Cl lg . b a m a x a0 ISTFI at ISTFI at ISTFI at ISTFI at Test a=b a=2b a=3b a=a0 1.00 2.4 2.4 7.39 1.88e2 7.46 2.29e2 C64 Table 3.10: Nonlinear (N = 3) streamfunction amplitude, first vertical mode. Parameter values: G = 500, e = 0.01, T f = 4 , and Cl lg . Chapter 3: Modified Rayleigh-Benard Convection ... 124 In order to determine the effect of the interaction of roll scales and imposed scales, numerical simulations based on the zeroth order equations (i.e. classical R-B convection) were conducted for G / G c = 1.004 (N44) and G / G c = 1.52 (N64) using identical initial conditions and parameter values as for C440 and C64 respectively. Results from N44 indicate that for the first vertical mode and a preferred wavenumber a p = 2.2, the ISTFI is 0.241, while for N64 the ISTFI at a p = 2.4 has a value of 245. Note (tables 3.7, 3.8, and 3.9) that for both values of the Grashof number linear and nonlinear results indicate that the ISTFI for the first vertical mode associated with the wavenumber a = a p decreases as b—>ap. Also note that although the nonlinear results for the first mode with G / G„ = 1.004 indicate that the wavenumber a = a n contains more energy than in the flat-plate case (N44), the opposite is true for numerical test conducted at G / G c = 1.52. However, an examination of mode 2 (results not presented) for G / G c = 1.004, at a = 2.2" indicates a ISTFI of ~ 0 for N44,.with 4.37e6 for C440 and ISTFI increasing as b -> a p . Similarly for G / G c = 1.52, at a = 2.4, the ISTFI is 3.03e"13 for N64 and 2.25e2 for C64. Thus the higher harmonics play an increasingly important role in the dynamics of the modified R-B convection model compared to the classical problem. This suggests the need to incorporate more vertical modes for numerical simulations based on the truncated equations as compared to the zeroth order equations (for the same value of the Grashof number). This is particularly true when resolving the thermal structure of the layer. The real and imaginary components of the streamfunction ( V | / a l , y a l ) and temperature (9 R , , 0 a l ) for the first vertical mode of C440 (with b = 1) are given in figure 3.16. From \ | / a l and 6 J , , the influence of T E R M (1) and T E R M (4) of (3.35) and (3.36) are readily identified (i.e. Chapter 3: Modified Rayleigh-Benard Convection . 125 0.02 0 -0.02 -0.04 -0.06 -0.08 -0.1 -0.12 Re(STF) Im(STF) 1 1 V 1 V 1 2 4 wavenumber CD -a 3 " a , £ cd 2 4 wavenumber Re(Temp) T Im(Temp) 2 4 wavenumber CD T3 3 £ cd -0.01 0 2 4 wavenumber Figure 3.16: Final spectra. Truncated equations test C440. Top left to right: Real and imaginary components of the streamfunction. Bottom left to right: Real and imaginary components of the temperature. Parameter values: G = 330 , N = 3, b = 1, e = 0.01, and T f = 12 with C l l f initial conditions. Chapter 3: Modified Rayleigh-Benard Convection 1 2 6 Xi B 3 C > 0.5 1 - 1 . 5 2 2.5 3 nondimensionalized time Figure 3.17: Linear and nonlinear wavenumber development. Truncated Equations. Wavenumber development based in linear test indicated by the solid line, with the dashed line representing nonlinear development. Note the change in scale associated with initial development at the imposed wavenumber and the shift in dominance to the harmonic at 2b. Also note that the nonlinear results lead the linear transition in scales. Parameter values: G = 500 , A = 10, b = 1, e = 0.01, and T f = 4 with CUg initial conditions. Chapter 3: Modified Rayleigh-Benard Convection ... 127 the delta function terms associated with a = b). Similarly, the shifting left and right in wavenumber space of the main lobe (TERM (2) of (3.35)) centered at 2.4 is apparent from all four figures in 3.15. Wavenumber development in physical space for the linear test C60 and the nonlinear test C64 with G / G c = 1.52, are presented in figure 3.17. As before, the flow is initially dominated by the forced scale associated with the imposed wavenumber 'b' (cf. figure 3.15). As the flow develops a transition to higher wavenumbers occurs (for this case) rather abruptly at x ~ 0.6. Thus as expected, when energy is made available to a wide range of wavenumbers, the transition to the preferred wavenumber occurs more rapidly (cf. transition at x ~ 1 for CIII initial conditions). The deviation of the nonlinear and linear cases for x > 3.4 is attributed to the nonlinear transfer of energy resulting from T E R M (3) of equation (3.35). Comparison of the results of modified and classical R-B convection Nonlinear simulations based on the truncated equations have provided information relating to the effects of a warped upper surface on scale selection. In particular, the results indicate a range of wavenumbers that is predicted to contain the dominant mode, and is centered around the preferred wavenumber associated with classical R-B convection. In order to further understand the effects of the warping on flow dynamics, the transport of heat by the roll circulation and the modification of the mean temperature profile were considered for both the classical and modified R-B domains for G / G c = 1.52 (test N64 and C64) with b = 1 and e = 0.01. (Note that values of the parameters for test N64 were the same as those for C64 which are given in table 3.4). Chapter 3: Modified Rayleigh-Benard Convection 128 The dimensional flux of heat from the surface to the fluid at a boundary may be calculated using icAT —~ 5Q = - ^ f [ ( V 0 ) n ] | z r l Z=0,(l+ecos(bX)) (3.55) where V 0 is the nondimensional temperature gradient, n is the unit vector directed into the fluid, the normal component of the temperature gradient is evaluated at the boundary, and the bar denotes averaging over the horizontal. Transforming (3.55) to mapped space, expanding to first order in e (the amplitude of the warping), and averaging over the entire horizontal domain, the heat flux into the fluid at the lower (8Q L ) surface, and out of the fluid at the upper surface (SQu) becomes 8 Q L = K A T 2H and a n k A T 2H 2-(l/2)Xk[e o R k -e0 b R k ] 2-(l/2)Xk(-l) k[e R k-ee R k] (3.56) (3.57) where 0 R k and 6 b k denote the value of the real component of the temperature in wavenumber space associated with a = 0 and a = b respectively. Based on (2.3) the average rate of change of temperature over the domain may be related to the heat flux at the boundaries by K A T (3.58) where A F = 2Fl[8Qu - 5 Q L ] / ( K A T ) is a nondimensional net heat flux across the layer. Thus the system will approach steady-state as [8Q U -8Q L ]—> 0, with the sign of (3.58) indicating Chapter 3: Modified Rayleigh-Benard Convection ... 129 whether or not the layer is heating (dT / dt > 0) or cooling (dT / dt < 0). Note that from (3.56) and (3.57), dQv = 8Q L when ^ , ( 9 R k - e 9 f k ) = 0 where the summation is only over the odd k=odd vertical modes. Finally, based on the convective time scale defined by (3.10) (i.e. [ H 2 / v ] ) and introducing a thermal time scale associated with (3.58) given by [2H 2 / ( K A F ) ] , the ratio of the convective to the thermal time scale becomes " 1 1 (3.59) ^ A F 2a Thus it will be assumed that the solution has reached steady state if (3.59) is less than some specified value, (say) 0.1. This means that the rate that the temperature of the layer (as a whole) is changing, is an order of magnitude less than the rate by which conduction is able to transport the heat through the domain. Example: Cases N64 and C64 Depicted in figure 3.18 are the ISTFI and the ITempI for G / G c = 1.52 based on N64 (top) and C64 (bottom) corresponding to the classical and modified domains respectively. Based on (3.4), AF ~ 0 for N64, and AF « 0.042 for C64, thus from (3.59) both of the test cases may be considered as having reached steady-state. Results indicate that for N64 the average heat flux through the domain is (16.2)KAT/ H , and (19.9)KAT for C64. Thus giving a relative increase of ~ 23% in the net heat flux through the domain at % - 4 (~ 2.9 hours) as a result of the warping of the upper boundary. Chapter 3: Modified Rayleigh-Benard Convection 130 250 2 4 6 8 wavenumber 2 4 6 wavenumber ISTFI 250 200 h 0 2 4 6 8 wavenumber -a 3 > s c3 50 40 h 30 h 20 10 h 0 2 4 6 8 wavenumber Figure 3.18: Final spectra for N64 and C64. Top left to right: ISTFI and ITempI associated with N64 (i.e. classical R-B convection). Bottom left to right: ISTFI and ITempI resulting from C64 (i.e. modified R-B convection). Parameter values: G = 500, N = 3, b = 1, and e = 0.01 with Cl lg initial conditions. Chapter 3: Modified Rayleigh-Benard Convection ... 131 Note, however, that the mean temperature profile illustrated in figure 3.19 shows little difference between the two cases presented. This discrepancy is attributed in part to the way in which the temperature profile and the heat flux values were calculated. Figure 3.19 has been constructed based on averaging the temperature field over the central eight rolls in physical space for both cases, while the heat flux calculations are based on the final wavenumber spectrum. Although the warping results in the simultaneous development of rolls throughout the entire domain at the imposed wavenumber b, the flow associated with the classical domain requires the propagation of the initial perturbation. Therefore, by the end of the model run there exists regionswithin the classical domain which remain undisturbed while there have been modifications to the temperature field throughout the entire domain associated with the corrugated boundary. 3.8.5 Summary and Discussion of Numerical Results The perturbation analysis of section 3.3.2 has resulted in an opportunity to study both the presumed effects of an induced gravity wave field, and the two-dimensional nonlinear exchange of energy alone as possible mechanisms for large aspect ration H R V . The latter is of interest in association with laboratory observations as the former mechanism is obviously not plausible in this case. In agreement with results of Getling (1983) and Chang and Shirer (1984), upscale energy transfer was observed only for a value of the Prandtl number of 0.1 with G / G c = 1.52. Results from test conducted within the laminar conductive regime, suggest an extension of the current analysis to three-dimensions in order to determine whether dimensional restriction of the flow accounts for discrepancy between two-dimensional results and experimentally observed scale selection above for G / G c > 1. Chapter 3: Modified Rayleigh-Benard Convection 132 nondimensional temperature Figure 3.19: Mean temperature profile. Classical versus modified R-B convection. The mean temperature profile for N64 (dashed line) and C64 (solid line) shows little difference between the two cases. See text for discussion of the heat flux. The short-dashed line represents the mean temperature profile associated with the conductive zeroth order state. Parameter values: G = 500 , N = 3, b = 1, and e = 0.01 with Cflg initial conditions. Chapter 3: Modified Rayleigh-Benard Convection ... 133 Combining the results of zeroth order and first order equations, the perturbation analysis suggests that an induced circulation will result for all nonzero values of the Grashof number at the imposed scale when G / G c < 1 and at multiple scales when G / G c > 1. Truncation of the perturbation analysis has lead to results which further suggest that competition to form the dominant scale exists between the wavenumber associated with 'free-flow' development (i.e. flow between parallel bounding surfaces) and the imposed scale. Dominance will be determined by the magnitude of the energy contained within each of these two wavenumber regimes. In agreement with results of the perturbation analysis, the truncated analysis indicates the existence of an induced circulation for all nonzero values of the Grashof number. However, results of the perturbation analysis have been shown to be misleading as the truncated analysis indicates that not only the imposed wavenumber but its harmonics, play crucial roles in induced H R V circulations. Within the conductive regime ( G / G c <1) , results indicate multiple scaling with the imposed wavenumber dominating. However in the convective regime ( G / G c > 1), the harmonic closest (in wavenumber space) to the preferred scale of free-flow development will contain the maximum amount of energy. Thus, dominance of a particular scale (i.e. forced or free) will be confined to a range of wavenumbers given by ( 2 / 3 < a / a p < 4 / 3 ) where a p is the preferred wavenumber associated with classical R-B convection. With respect to heat flux estimates, results for G / G c = 1.52 indicate a 23% increase in the average heat flux through the modified domain, compared to average heat flux through the classical domain. These results are attributed to an increased importance of vertical modes greater than one when the domain is modified to include a corrugated upper surface. Chapter 3: Modified Rayleigh-Benard Convection ... 134 3.9 Application of Results to the Atmosphere Numerical simulations presented in section 3.8 focused on values of the Grashof number for which laboratory experiments have observed two-dimensional steady-state convection (see table 2.1 where Ra = 2 G G ) . Results (section 3.8.4) indicate that the dominant scale is limited to a band of wavenumbers (namely 2 / 3 < a / a p < 4 / 3 , where ap is the preferred wavenumber associated with classical R-B convection). More specifically, as the flow develops, the harmonic of the imposed wavenumber 'b' closest to ap (in wavenumber space) dominates. As the present study was motivated by observations of H R V within the PBL, the current section presents an investigation into roll formation at geophysically observed values of the Grashof number (i.e. 5000 to 25,000, based on the current definition of the Grashof number (3.11), an eddy viscosity and eddy conductivity of -50 m2/s giving an eddy Prandtl number of 1, and observations presented by Miura (1986)). Laboratory observations presented in table 2.1 indicate that this range of Grashof numbers lie within the turbulent convection regime. Thus it is anticipated that flow development will be associated with roll formation which is both time-dependent and aperiodic. 3.9.1 An Atmospheric Test Case In order to investigate the influence of a gravity wave field on boundary layer dynamics, numerical simulations were conducted for flow associated with both the classical (isolated boundary layer) and the modified (coupled boundary layer-gravity wave) domains. Al l tests were conducted using G / G c = 22.8 (i.e. G = 7500), an eddy Prandtl number (cr E) of one, and a final nondimensional time of 0.5. Spectra were initialized using class II initial Chapter 3: Modified Rayleigh-Benard Convection ... 135 conditions with S = 0.001, s = 0.5 and a = 2.2. (The initial spectrum for the ISTFI is given in figure 3.20). Thus the initial perturbations were weak, and centered about 2.2 in wavenumber space (i.e. a c ) . A (relatively) broad initial bandwidth was prescribed so that energy would be available to a wide range of scales. Figure 3.20: Initial spectra. Magnitude of the streamfunction. Classical R-B convection Numerical simulations based on the zeroth order equations (which govern the dynamics of classical R-B convection) were a logical first step since the preferred scale of convection would provide information used to estimate the response of the flow associated with the modified domain, to the choice of imposed wavenumber 'b'. Also, results from the zeroth order equations may be compared with those from other 'isolated' boundary layer studies at high Grashof Chapter 3: Modified Rayleigh-Benard Convection ... 136 numbers (e.g. Rothermal and Agee, 1986; Sykes and Henn, 1988). Model parameters used in the numerical simulations are presented in table 3.11. Test a G Tf N A nodes step size method A4 1 7500 0.5 3 15 151 5xl0" 4 rk4. CII A l 1 7500 0.5 5 15 251 lxlO" 3 rk4, CII A7 1 7500 0.5 7 15 151 5xl0" 4 rk4, CII Table 3.11: Zeroth order equations. Numerical test parameters. Illustrated in figure 3.21 is the final ISTFI spectra of A7 with 7 vertical modes. For the duration of the simulations, the dominant range of scales was found to be relatively insensitive to the number of vertical modes, with test cases A4, A l and A3 all indicating a preferred wavenumber of ~3. This suggests that scale selection is adequately modeled and that roll formation is dominated by horizontal scales which are smaller than those observed at the critical Grashof number (i.e. to the left of a c (-2.221) in wavenumber space). The streamfunction in physical space for A7 is also given in figure 3.21. However, unlike scale selection, sensitivity to the number of vertical modes is clearly evident in the resolution of the thermal features (figure 3.22). The top figure is for N = 3, the lower is for N = 7. In particular, note the evidence of updrafts and downdrafts within the flow which are characterized by deformation of the isotherms. In figure 3.22 (lower) the deformation of the isotherms extends from one surface (almost) to the other resulting in a compacting of the isotherms in a thin region near the boundaries, known as the thermal boundary layer. Chapter 3: Modified Rayleigh-Benard Convection . 137 ISTFI 600 -500 -co 1 400 -f* 300 -0 2 4 6 8 10 12 wavenumber X -» Figure 3.21: Atmospheric test case. Classical R-B convection. Top: ISTFI. Bottom: Streamfunction in physical space. Note that the preferred scale of convection lies in the range of [2,4] and is approximately 3. Parameter values: G = 7500, G E = 1, N = 7 with CH initial conditions. Chapter 3: Modified Rayleigh-Benard Convection . 138 -10 -5 0 5 10 X -> -10 -5 0 5 10 X -» Figure 3.22: Temperature in physical space. Top: Temperature in physical space for A4 with 3 vertical modes. Bottom: Temperature in physical space for A7 with 7 vertical modes. Note the regions of nonphysical temperatures (i.e. areas where ITI> 1 as denoted by the arrow in the top figure) resulting from inadequate vertical resolution. Warm updrafts and cool downdrafts are evident in the lower figure (the arrow denotes a warm updraft). Also note compacting of the isotherms near the boundary. (See text for details.) Chapter 3: Modified Rayleigh-Benard Convection ... 139 At subcritical values of the Grashof number, the isotherms are parallel to the boundaries and therefore the depth of the thermal boundary layer is equal to the distance separating the two plates. However, as the Grashof number increases above the critical Grashof number, the roll circulation results in a deformation of the isotherms and a positive vertical heat flux, thus forming warm updrafts and cold downdrafts. Since the magnitude of the circulation increases as G increases, the thickness of the thermal boundary layer will decrease. If the number of vertical modes is insufficient, nonphysical regions (typically) within the thermal boundary layer will develop where the temperature is greater (or smaller) than the surface temperature (note the upper figure of 3.22). Therefore, the existence of these nonphysical regions near the surface (rather than the range of dominant wavenumbers) provides a good guideline by which to determine whether or not the numerical solution obtained accurately represents the solution of the physical model. This naturally suggests that an increasing number of vertical modes will be required as the Grashof number increases. (Compare the upper and lower figures of 3.22). In particular the well-mixed PBL which is characterized by a constant mean (potential) temperature profile (see for example Stull, 1988) above the surface layer, is analogous in the current model to a mean temperature profile which is zero throughout the domain except in the thermal boundary layers where the temperature will increase (decrease) to satisfy the lower (upper) boundary conditions. The mean temperature profile for test A7 is given in figure 3.23. Note the formation of an isothermal middle layer and an upper and lower thermal boundary layers. This suggests, that the model is able to produce results which compare well with observations of mean convective boundary layer characteristics. Chapter 3: Modified Rayleigh-Benard Convection 140 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 nondimensionalized temperature Figure 3.23: Mean temperature profile. Classical R-B convection. The mean temperature profile (solid line) was calculated by averaging across the middle 8 rolls depicted in figure 3.21. Note that the mid-section of the region between the plates will tend to a vertical linear profile (i.e. a well-mixed layer) as the number of vertical modes included in the numerical model increases. Also note the thermal boundary layer near the upper (rj = 1) and lower (t| = 0) surfaces. As the thermal boundary layer thickness decreases within increasing Grashof number, the number of vertical modes required to resolve the thermal boundary layer will also increase. Note that the dashed line represents the linear mean temperature profile associated with the conductive state. Chapter 3: Modified Rayleigh-Benard Convection ... 141 (As noted in section 2.5.1, the assumption of a perfectly conducting, isothermal upper boundary condition is highly idealized and results in an thermal boundary layer at the upper plate which is not observed in the real atmosphere.) Modified R-B convection Dynamics of the coupled system were modeled using the truncated equations. Based on the results of the isolated boundary layer, it was assumed that a p e [2,4] with ap estimated to be approximately 3. Results of chapter 4 (which were based on numerical simulations of Clark et al. (1986)) suggest an imposed wavenumber of b = 2.1. However, since the imposed and preferred wavenumbers are not significantly different, numerical test were also conducted for b = 1.6 and b = 1. Note that an amplitude of the warping of 1% of the boundary layer depth (i.e. 8 = 0.01) was assumed. Remaining parameter values are given in table 3.12. Test c G b T f N A nodes step size method A3 1 7500 2.1 0.5 3 15 151 5xl0" 4 rk4. CII A30 T 7500 2.1 0.5 5 15 151 5xl0 ' 4 rk4, CII A2 1 7500 1.6 0.5 3 15 151 5xl0" 4 rk4, C H A5 1 7500 1.0 0.5 3 15 151 5xl0" 4 rk4, CII A6 1 7500 1.0 0.5 5 15 151 5xl0" 4 rk4, CII A8 1 7500 1.0 0.5 7 15 151 5xl0" 4 rk4, CII Table 3.12: Truncated equations. Numerical test parameters. If we assume that the results from the truncated equations for the laminar convective regime apply to all convective regimes, then it is expected the effect of the warping on flow dynamics will be to isolate (within a range of wavenumbers) the preferred scale associated with classical R-B for the Chapter 3: Modified Rayleigh-Benard Convection ... 142 same value of the Grashof number (assuming that a such a preferred range of scales exits). Thus it was anticipated that for b = 2.1 the forced wavenumber would dominate or possibly the first harmonic a = 2b = 4.2; for b = 1.6 the first harmonic is expected to dominate i.e. a = 2b = 3.2; and finally, for b = 1 a dominant harmonic of 2, 3, or (possibly) 4. These results would help to identify the preferred scale(s) from the broad spectrum of wavenumbers [2,4] containing significant amounts energy in the flat-plate case. Also, based on results from the classical model it was not anticipated that vertical mode resolution would bias horizontal wavenumber selection. Comparison of figures 3.21 (classical R-B convection) and 3.24 (modified R-B convection) indicates that for b = 2.1 the dominant band of wavenumbers has shifted slightly to smaller values (i.e. larger scales) and is now approximately centered around the imposed wavenumber. Results for b = 1.0 and b = 1.6 (figure 3.24) were not as expected, based on the assumption that the preferred wavenumber was ~ 3. The middle figure of 3.24 corresponding to an imposed wavenumber b = 1.0 clearly identifies the first harmonic (i.e. a = 2) as dominant. These results are somewhat surprising since figure 3.21 suggests that the third harmonic would play a more significant role. Results from b = 1.6 also suggest that the preferred scale is associated with larger scales than expected. Combining the range of scales over which these two case would dominant, the overlapping range of wavenumbers suggests that the preferred scale lies in the range [1.5,2.4]. Also note that the ISTFI associated with b = 1.6 (A2) is approximately twice the value of the ISTFI for A30 with b = 2.1. This suggests that in fact, the preferred wavenumber may be closer to 1.6 than 2.1. However, further testing would be required in order to determine whether or not this is in fact the case. Chapter 3: Modified Rayleigh-Benard Convection . 143 CO T3 3 3 CD T3 3 4 6 8 10 12 14 wavenumber Figure 3.24: Final spectra. Magnitude of the Streamfunction. Top: b = 2.1 with N = 5 . Middle: b = 1.0 with N = 7 . Bottom: b = 1.6 with N = 3. Chapter 3: Modified Rayleigh-Benard Convection ... 144 Results from the truncated equations raises questions as to the ability of the two-dimensional nonlinear interactions to produce scales which are larger than Xc, particularly at high values of the Grashof number. Although the present results do not contradict results from test conducted in the laminar regime (section 3.8.4), they do suggest the need for further study into nonlinear scale selection associated with classical R-B convection. Recall (section 2.2) that Rothermal and Agee (1986) found an upscale energy transfer for large values of the Grashof number. However, Sykes and Henn (1988) found no such increase in scale. Further study is warranted. In summary, model results based on the truncated equations have indicated roll formation at wavenumbers (1.6) which are less than the critical wavenumber (2.221), and correspond to an aspect ratio of ~ 4 in physical space. However, forcing associated with b = 1 did not result in dominant large scales corresponding to an aspect ratio of ~ 6. Instead, aspect ratio rolls of ~ 3 associated with the first harmonic of the imposed wavenumber were found to dominate. 3.9.2 Discussion of Results Based on observations of Miura (1986), an eddy Prandtl number of 1 and the definition of the Grashof number, the range of Grashof numbers for which there have been observations of H R V was estimated to be [5000,25000]. Numerical test were conducted using G = 7500. In order to determine the effect of the warped upper surface on flow dynamics, two types of numerical simulations were conducted. The first were based on the zeroth order equations (i.e. classical R-B convection) while the second type of numerical simulations were based on the truncated equations (i.e. modified R-B convection). Chapter 3: Modified Rayleigh-Benard Convection ... 145 Results of the flat-plate (or isolated boundary layer) case, indicate that the model is able to reproduce physically realistic results with respect to the mean thermal properties of the convectively driven PBL. In particular, a well-mixed layer occupying ~ 70% of the domain and a thermal boundary layer (or surface layer) near the bounding surfaces. Al l numerical tests conducted based on the truncated equations, assumed an amplitude of the warping of the upper surface of 0.01. The choice of 8 was arbitrary and was felt to be an 'order of magnitude' estimate of the potential of the upper flow to warp the mean temperature field. Numerical modeling results presented by Carruthers and Moeng (1987) suggest an 8 of approximately 0.025 (see figure 1 of Carruthers and Moeng, 1987). Although numerical simulations presented from the current model are for e = 0.01, a test conducted with G = 7500, b = 1, o E = l , N = 5, A = 15 , T f = 0.5 and 8 = 0.05 (compare with test A6, and the middle figure of figure 3.24) indicated the dominance of the second harmonic in agreement with results conducted with e = 0.01. Results from the truncated equations suggest that the preferred scale (of the flat-plate case) lies in the wavenumber range [1.5,2.4]. This result was unexpected, since the flat-plate case indicated a preferred scale closer to 3. Further investigation is necessary in order to determine whether the transition to smaller wavenumbers would result if the flat-plate case is run for a longer period of time. Consideration of the net average heat flux transported by the roll circulation, resulted in a 27% increase associated with the warped circulation (test A8 with a mean heat flux based on (3.56) and (3.57) of ( (61.21)KAT / H ) , compared to the flat-plate case (test case A7 with a mean heat flux of (48.25)KAT / H) . Paralleling laminar results (section 3.8.4), this increase in efficiency Chapter 3: Modified Rayleigh-Benard Convection ... 146 is attributed to an increase in the role of the convective scales, particularly in association with the higher vertical modes. The question of the partitioning of the energy between the various terms of the kinetic energy budget is currently unresolved. However, such an analysis is of interest as quantification of the value of each term of the kinetic energy budget, would help to better understand the effect of the warped upper surface on the energetics of the system. Comparison of the thermal and convective time scales based on (3.59) indicate that the flow associated with both the zeroth order equations (AF ~ 1.92) and the truncated equations (AF = 0.46) have not yet reached steady-state. However based on calculations of the net heat flux at the final value of the nondimensional time, the modified system appears to be closer to steady-state than the classical R-B system. This agrees with the suggestion (above) that the preferred wavenumber associated with the flat-plate case may not have reached its final value. Chapter 4 The Stable Upper Layer and Induced Gravity Waves As discussed in section 2.5.3, the effects of boundary layer dynamics on the overlying stable layer are incorporated by assuming that the effect of the periodic initial warping of the inversion by the boundary layer roll vortices is analogous to flow over topography. The characteristics of the stable upper layer are investigated using stability and mean velocity profiles estimated from temperature and velocity profiles of Clark et al. (1986) which they used in their two-dimensional, finite-difference model. The profiles are presented in section 4.2. The stable upper layer is modeled using both a one-layer (section 4.4), and a threQ-layer (section 4.5) model. Linear perturbation analysis (section 2.1) is used to identify the resonant mode(s) of the atmosphere as well as the principal scale of influence at the inversion. 4.1 Description of the Physical Model Consider a scenario as depicted in figure 4.1 (cf. figure 2.3). The convective activity has just filled the domain of the boundary layer and the warping of the inversion is proposed to occur at the preferred scale of the roll vortices (as modeled by the classical R-B convection), and is given by X = 2K I a 0 . It is suggested that the initial value of the wavenumber will be associated with the critical wavenumber (i.e. a 0 = a c ) where a c has a nondimensional value of 2.221 (section 3.6.1). If the inversion is warped then the mean wind shear will result in the generation of gravity waves in the stable upper layer at the forced scale. 147 Chapter 4: The Stable Upper Layer and . 148 U 27t/a 0 < > z=0 inversion CBL Figure 4.1: Topography of lower surface of the stable layer. Initial warping of the inversion by the convectively driven roll vortices will occur at the critical wavenumber ao. Note that the vertical arrows denote updrafts and downdrafts associated with the roll circulation. The inversion, assumed rigid, is represented by H(x) = h 0 sin(a0x) where h 0 is the amplitude of the warping and is assumed known. Numerical results of section (3.8) indicate that as a result of the nonlinear transfer of energy, there exists energy (within the roll circulation) not only at a single scale, but in a broad spectrum of scales. Thus, although the warping of the inversion by the convective activity is modeled for simplicity at a single length scale, it should be noted that the warping actually occurs on a large number of scales with various amplitudes. Thus the 'topography' should be assumed of the form However, since the response of the upper layer is investigated using a linear analysis, the superposition principle allows the investigation of each of the individual wavenumbers in order to develop the net response. 0 Chapter 4: The Stable Upper Layer and ... 149 Of particular interest is the resonant mode which is predicted (by linear theory) to have a steady-state amplitude which is unbounded. Although the spectral energy associated with the roll circulation will be banded about the preferred scale, if sufficient energy (in the roll circulation) exists at the wavenumber corresponding to the resonant mode of the upper layer, it is plausible that the resonant mode will develop with the maximum growth rate. It is then further hypothesized that the maximum rate of nonlinear transfer of energy to other gravity wave scales will occur at the height corresponding to the height of maximum energy associated with the resonant mode. Finally, in order to determine the mode of maximum influence at the inversion, the transfer of energy from the resonant mode to other scales is assumed to follow an equal partitioning of kinetic energy. Thus the general (linear) solution is scaled based on kinetic energy assumptions, and the mode with maximum vertical amplitude at the inversion is assumed to have maximum influence on the dynamics of the convectively driven planetary boundary layer. (Note that such a scaling is highly idealized but necessary in the absence of a nonlinear analysis.) 4.2 Atmospheric Temperature and Velocity Profiles of Clark et al. 1986. The choice of vertical profiles used in the numerical simulations of Clark et al. (1986) (referred to from now on as C86), was motivated by observations of Kuettner et al. (1987) (or K87) as part of the 'Convection Wave Project' launched by the National Center for Atmospheric Research (NCAR) in Boulder, Colorado. The observational study was designed to investigate the nature of the gravity waves which are generated by convective activity (as opposed to 'mountain waves' Chapter 4: The Stable Upper Layer and ... 150 formed by topography) using instrumented aircraft, and of particular interest are the observations from the 12 th of June 1984, which were obtained over the western plains of Nebraska. The velocity and temperature profiles used by C86 are presented in figure 4.2. Note that the mean wind shear is assumed constant with height with a value of 7xl0~ 3 s _ 1 over the first 3 km, with a constant mean wind of 10 m/s over the remainder of the domain. From the temperature profile, it can be seen that the numerical model extends into the lower portion of the stratosphere which was considered important as it was felt that the tropopause may provide an important source of internal wave reflection, and that domain size would influence scale selection. Using a horizontal domain of 30 km and a vertical domain of 15 km extending into the stratosphere (0-2 km boundary layer, from 2-15 km as the upper layer), numerical results of C86 indicate an average gravity wavelength of 11 km just below the tropopause, and approximately 7.5 km in the lower part of stable layer near the inversion. Field observations from K87 taken at a height of approximately 1 km above the inversion indicate an average wavelength of 6 km. In order to develop a theoretical model of the atmosphere, the environmental temperature profile and the form of the velocity profile assumed in the two-dimensional numerical test of C86 were used to estimate potential temperature and atmospheric stability profiles which are presented in figure 4.3. Aside from the 'peak' in the square of the Brunt-Vaisala frequency (N 2) profile at the height of the boundary layer (~2 km), the stability may be seen to decrease gradually over the troposphere, with the tropopause clearly evident at a height of 9 km. Note that within the C B L (indicated by the line at 2 km which refers to the height of the convective boundary layer up to the capping inversion) only a portion of the N 2 profile (based on potential temperature) indicates instability. In the presence of cumuli (which was the case for the temperature profile of figure 4.2) er 4: The Stable Upper Layer and 16 oo X B J4 .£? '5 X Figure 4.2: Velocity and temperature profiles from Clark et al. (1986). 16 | 1 1 1 16 14 - / - 14 12 \- / - | 12 10 h 10 h op 250 300 350 degrees (K) 400 0 -0.0001 0.0002 (sec)A(-2) 0.0005 Figure 4.3: Estimated profiles of potential temperature and square of the Brunt-Vaisala frequency. Chapter 4: The Stable Upper Layer and ... 152 the potential temperature is not constant with height within the convectively driven P B L but increases slightly in the upper part, thus K87 suggests the use of a near constant equivalent potential temperature as a more applicable definition of the C B L rather than one based on potential temperature. The wavelength associated with the Scorer parameter (defined below (4.2), and illustrated in figure 4.4) gives an estimate as to the extent of vertical propagation of waves of horizontal wavenumber k based on the stability and mean wind profiles of figure 4.2. Comparison of numerically obtained results of C86 and theoretically predicted values of the resonant and dominant modes from a one-layer model of the atmosphere is presented in section 4.4, while the resonant mode of a three-layer model is presented in section 4.5. 16 14 h 12 10 0 0 5 10 15 20 (km) Figure 4.4: Wavelength associated with the Scorer parameter. Chapter 4: The Stable Upper Layer and ... 153 4.3 Governing Equations Consider two-dimensional, inviscid flow and assume that the flow variables may be separated under Reynolds decomposition into mean and perturbation components of the form f(x,z) = F(z) + f (x,z) where f represents the velocity components (u,w), the dynamic pressure p, the density p, or the potential temperature 9. (Note that the potential temperature is defined by R / c 9 = T(p s / p) p where p s is a reference pressure, T is the temperature, R is the ideal gas constant and c p is the specific heat at constant temperature. See for example Holton, 1992.) In particular, assume p(t, x, z) = p 0 + p' (t, x,z) p(t,x,z) = P(z) + p' (t,x,z), 9(t, x, z) = 9(z) + 9' (t, x, z), u(t,x,z) = U(z) + u'(t,x,z), and w(t,x,z) = w'(t,x,z). (4.1) Then the equations governing the perturbation quantities become (e.g. Holton, 1992) 3u' -au' , 9U 1 dp' „ — + U — + w' — + — - ^ = 0 dt dx dz p o dx 3w' -3w' 1 dp' 0' —— + U - — + f--=g = 0 dt dx p0 dz 9 du' aw' + = 0 dx dz 99' -39' ,39 n — + U — + w' — = 0 (4.2) 3t 3x 3z 9' D' where we have used — ~ ——. If it is further assumed the flow is time-independent (i.e. fi Po stationary with respect to the ground), the equation governing the vertical perturbation velocity w' (x, z) may be written Chapter 4: The Stable Upper Layer and . 154 av a 2 w + l 2 w' = 0 (4.3) dz2 where 1 is the Scorer parameter (units rn1) defined by l 2 = N2 (z ) 1 d 2 U(z) U 2 (z ) U(z) dz 2 and N 2 s g[d(ln 0) / dz] is the Brunt- Vaisala frequency. 4.4 A One-Layer Model of the Stable Upper Layer. With the height coordinate z having its origin at the top of the boundary layer, the one-layer model of the atmosphere assumes a constant N 2 and mean density1 with height, and a linearly increasing mean velocity profile given by U(z) = U( l + z / L ) , as illustrated in figure 4.5. This constant shear profile is assumed to extend throughout the model domain defined by z e [0,°°) and at not just from the inversion to a depth of 1 km, as was used by C86 (c.f. figure (4.1)). Assuming a form for the vertical component of the perturbation velocity, w' (x,z), given by w' (x,z) = w(z)sin(kx), (4.3) may be simplified to give the following equation governing the vertical dependence w(z) where we have introduced Z = (l + z / L ) , a 2 = k 2 L 2 , and Ri = r— as the mean Richardson U 2 number. Letting |1 = ^ R i - 1 / 4 , the general solution bounded as (ocZ) —> °° is given by 1 In theory, a constant value of N 2 implies a vertical density profile which decreases exponentially with height at a rate of (-N2/g). However for the one-layer and three-layer models used in the current study, it is assumed that the variations of density with height are negligible, after the Boussinesq approximation. (See appendix I, case 4, for a discussion on the inclusion of density variations after Wurtele et al. (1987).) N 2 L 2 Chapter 4: The Stable Upper Layer and . 155 Figure 4.5: One-layer model of the stable layer. Assumed vertical profiles of mean wind (U), the square of the Brunt-Vaisala frequency ( N 2 ) , and potential temperature (0) profiles for the one-layer model. Note that the mean wind increases linearly with height throughout the entire domain defined by z e [0, °°). Also N 2 is assumed constant with height and since N 2 = g(d ln(0) / dz), the vertical potential temperature profile increases exponentially with height given by 0 = 0O exp(N 2z / g). Letting \i = V R i - 1 / 4 , the general solution bounded as (aZ) —> oo is given by w(Z) = A Z 1 / 2 K i f l ( a Z ) where ' A ' is an arbitrary constant and K^(ocZ) is related to the Bessel function of the First Kind of imaginary order and imaginary argument and is defined by (see Abramowitz and Stegun (1964) and Morgan (1947), v , i7TMx ) - I_ i ( i (x ) ^ 2 sinh(|i7X) or in integral form as K i t l (x ) = Jcos(ii.t)exp(-xcosh(t))dt. Thus, the corresponding general form 0 of the vertical perturbation velocity is given by Chapter 4: The Stable Upper Layer and ... 156 w' (x, Z) = A Z 1 / 2 K i u (aZ) sin(kx). (4.4) The horizontal component of the perturbation velocity may be found using the equation of continuity dw' (x, z) ro  ,  u (x,Z) = - dx J A T dz giving 1 a[z , / 2K i ( 1(aZ)l u'(x,Z) = - A - t ^cos(kx). (4.5) a dZ Finally, the particular solution may be found by applying the lower boundary condition which for topography of the form F£(x) = h 0 sin(k0x) may be written w (x,0) = Dt - u l -o = _ u h ° k o s i n ( k ° x ) • ( 4 - 6 ) z=o z dx Solving for the coefficient A using (4.4) and (4.5), the particular solution is given by w' (x,Z) = - J J k ° h ° Z 1 / 2 K i / a 0 Z ) s i n ( k 0 x ) (4.7) where a 2 = k 2 L 2 . From (4.7), it may be seen that given a warping of the inversion at a wavenumber k0, such that a o is a root of K; , i.e. K i t t ( a o ) = 0, the particular steady-state solution to the vertical component of the perturbation velocity predicts an unbounded amplitude for all values of Z greater than one. This value of the wavenumber(s) is the resonant mode(s) of the system. Chapter 4: The Stable Upper Layer and ... 157 4.4.1 Resonant Mode For values of the parameters estimated using the environmental profiles of C86, presented in section 4.2, the value of the resonant wavenumber were found by identifying the root of K ; (oc) = 0. The values of the parameters were estimated from values just above the 'spike' shown in figure 4.3 depicting the N 2 profile. Thus the mean wind at z = 0, i.e. U of = 5m / s , the vertical length scale L is estimated to lie in the range of [0.82,1] kilometers, and the value of N 2 is taken as l.lSxlO^s" 2 . Thus the Richardson number is estimated at have a value between [3.09,4.6], and the corresponding range of [L is given by [1.69,2.09]. Using M A P L E , with u. = [1.69,2.09] the only root of K i ( l(cc) found between [-0.3,14] was a r e s =[0.2995,0.4919] (where the subscript 'res' denotes resonant value). For an L of [0.82,1], these values of a are equivalent to a resonant wavelength between [17.20,12.77] kilometers. Although (4.7) predicts an infinite value for the vertical perturbation velocity associated with the resonant mode, it should be noted that the temporal development of the physical system will not realize such a solution as the nonlinearity of the system will result in a transfer of energy and the ultimate dissipation of energy at smaller scales. It is however, anticipated that the largest growth rate will be associated with the resonant mode as predicted by linear theory. However, the current analysis focuses on the behavior of the system for values of time for which linear analysis is assumed applicable, i.e. for t = e « 1. The vertical amplitude profile of w(Z), u(Z), and the kinetic energy profile (all of which have been divided by the coefficient A), corresponding to the resonant mode with Richardson number of 4.25, an L of 1 km, \i = 2, a r e s = 0.448, and resonant wavelength of 14 km, are illustrated in figure 4.6. The difference in the magnitude of the two perturbation velocity components may be seen from figure (4.6) where u m a x (Z) = 2.6wm a x (Z). Chapter 4: The Stable Upper Layer and . 158 -0.1 0 0.1 0.2 0.3 0.4 0 0.05 0.1 0.15 nondimensionalized amplitude nondimensionalized amplitude Figure 4.6: Resonant mode. One-layer model of the stable upper layer. Left: Profile of the vertical functional dependence associated with the perturbation velocity components w ( Z ) / A (solid line) and u ( Z ) / A (dashed line), where w(Z) / A = Z 1 / 2 K i M ( a r e s Z ) , u ( Z ) / A = (l /a)[d(Z 1 / 2 K^(a r e s Z)) /3z], and the subscript 'res' denotes evaluation at the resonant wavenumber. Right: Vertical profile of the horizontally averaged kinetic energy (4.8). Parameter values: L = 1, Ri = 4.25, p, = 2, and a r e s = 0.448075 . Chapter 4: The Stable Upper Layer and ... 159 This dominance in magnitude of the horizontal component results in a kinetic energy profile which is maximum at the inversion, where u(Z) attains it's maximum value, while the vertical velocity component reaches a maximum at a height of approximately 2.9 kilometers. Also note, that the vertical perturbation velocity associated with the resonant mode has a zero amplitude at the interface where Z = 1. Therefore, although it is expected that the resonant mode will be the primary mode observed throughout the domain, it is not expected to alter the length scale at the bounding surface. 4.4.2 Dominant Mode: Kinetic Energy Considerations From figure 4.6, we can see that the resonant mode is associated with a value of the vertical perturbation velocity at the inversion (Z = l) given by w(Z) = 0. Therefore, although the resonant mode will contain the maximum amount of energy, it does not have the potential to warp the inversion. Identifying the mode with the maximum warping potential, requires determining the value of the amplitude coefficient A of (4.4). Since the resonant mode contains a large amount of energy, it is assumed that this energy will be made available to the full spectrum of wave scales via nonlinear interactions. The choice mode here is to normalize the amplitude coefficient such that different gravity wave scales contain equal amounts of kinetic energy. Such an approximation is applied as in the absence of a fully nonlinear model of the stable upper layer there is no way to calculate the coefficient. We shall see (figure 4.8) that the dominant scale is well defined and is thus relatively insensitive to the choice of scaling. The expression for the horizontal component of the perturbation velocity (4.5), may be written explicitly as Chapter 4: The Stable Upper Layer and 160 u' (x,Z): kL (1 / 2 )Z~ 1 / 2 K U l (oZ) + Z 1 / 2 — ^ — -oZ cos(kx). Then the horizontally averaged kinetic energy may be calculated using 2K/k < K E > = — f((u') 2 +(w') 2)dx 4TC J (4.8) to give < K E > = ^ [ F a ( Z ) ] . (4.9) where F a (Z) = [z , / 2K i/aZ)] 2+-J 2- r 1/2 dK i | t (ctZ)" ( l / 2 ) Z - ' " K 1 ( l ( a Z ) + Z  ^ (4.10) From figure 4.6 it can be seen that the kinetic energy associated with the resonant mode is maximum at the inversion. Therefore, in order to scale waves of the form (4.4) it is assumed that the kinetic energy is partitioned equally between all wavenumbers at the inversion where Z = 1. Equation (4.9) may then be scaled to (say) 1/4 (eliminating the factor of 4) and the coefficient A may be written 1 A = [F«(D] 1/2 Thus the expression governing the vertical component of the perturbation velocity (4.4) becomes w ' s c e a (x,Z) = | Z 1 / 2K i t l(ocZ)sin(kx) (4.11) where the subscript denotes the 'scaled' solution of the vertical perturbation velocity. Chapter 4: The Stable Upper Layer and ... 161 In order to identify the dominant mode at the inversion, consideration must be given to the relationship between the vertical perturbation velocity and the (potential) amplitude of the resultant warping of the inversion. Determining the mode which maximizes the height perturbation of the inversion will fix both the dominant scale and it's amplitude (i.e. the imposed wavenumber 'b' and warping of the inversion 'e', of chapter 3.) Analogous to (4.6) the displacement (rj) of the inversion from it's initial, undisturbed height at z = 0 may be written W(x,0) = ^ 2> (4.12) and since gives dt dri(x,0) 9ri(x,0) ^ ^ ^ , 0 ) dt dt 3x ri(x,0) = ( l / U ) J w ' (x,0)dx or r|(x,0) « (w' (x,0) / kU). (4.13) Of particular interest is the wavenumber associated with the maximum amplitude of warping at the inversion (i.e. the maximum value of T|); equation (4.10) gives r|(oc) «= w a ( l ) / k where the subscript on w emphasizes dependence on alpha. Therefore, identifying the dominant mode at the inversion is equivalent to identifying the wavenumber for which the expression (1 / a)[K^ (a) / [F a (1)]1/2 ] (or equivalently w' (x,0) / k) is a maximum. Again using M A P L E with [I = [1.69,2.09], only one maximum was found in the range of a from [0.3, 8] corresponding to Chapter 4: The Stable Upper Layer and ... 162 a d o m = [0.872,1.157]. For values of L between [0.82,1.0] km, this value of a is equivalent to a dominant wavelength between [5.91,5.43] kilometers. The vertical functional dependence for the perturbation velocity components of the dominant mode (with L = 1, Ri = 4.25, \i = 2, ocd o m = 1.094 and corresponding wavelength of 5.74 km) as well as the vertical dependence of the kinetic energy are plotted in figure 4.7, and note that the kinetic energy attains its maximum value at a height of approximately 0.1 km above the inversion. Note that the dominant mode is trapped at lower levels in the atmosphere when compared to the longer wavelength resonant mode depicted in figure 4.6. Figure 4.8 depicts the value of the w ( Z ) / a as a function of a (or equivalently as a function of the wavenumber) evaluated at the inversion. Recall that T|(a) °c w a (1) / k represents the deflection of the inversion by a wave of length 27t/k. Figure 4.8 also indicates that the wavelength associated with the dominant mode is well defined. Note: The normalization of the modes given by (4.4) could have been done using the kinetic energy contained by the particular solution (4.7) where it is assumed that the partitioning of energy to is such that all wavenumbers contains the same fraction (p) of the energy as contained in the forced mode at the inversion, i.e. < K E >= (3 < K E > p a r t i c u i a r . Then from (4.6) and (4.9) [ F a o ( l ) ] a n d < K E > = ^ [ F a ( l ) ] < K E > p a r t i c u l a r " K f c « x 0 ) Chapter 4: The Stable Upper Layer and . 163 N 5 1 1 / 1 1 1 t — / / / / / / — / / / / / 1 \ ~~~~~l -0.05 0 0.05 0.1 nondimensionalized amplitude 0 0.0025 0.005 0.0075 0.01 nondimensionalized amplitude Figure 4.7: Dominant mode. One-layer model of the stable upper layer. Right: Profile of the vertical functional dependence associated with the perturbation velocity components for the dominant mode. Solid line: Vertical perturbation velocity component. Dashed line: Horizontal perturbation velocity component. Left: Vertical profile of the horizontally average kinetic energy (4.8). Comments as in figure 4.6. Parameter values: L = 1, Ri = 4.25, \i = 2, and Chapter 4: The Stable Upper Layer and . 164 XI =s "o. E CO CD N c o CD E x> c o Figure 4.8: Values of T]^^ = [w a(l)] / a as a function of a. Plotted is the value of the horizontally averaged vertical perturbation velocity component, which has been scaled using kinetic energy considerations, as a function of a. (Equation (4.11)). Note that the maximum deflection at the inversion is associated with the dominant mode at ocd o m = 1.094. Parameter values as in figure 4.8. gives a value of the amplitude coefficient A P ' / 2 U k 0 h 0 K,(a 0) 1/2 (4.14) and finally an expression for the vertical component of the perturbation velocity as w' (x,Z) = P "2 U k 0 h 0 K , ( a 0 ) F a 0 (D F a(D 1/2 [z , / 2 K i M XaZ)]sin(kx) (4.15) Chapter 4: The Stable Upper Layer and 165 Using the notation introduced by (4.11), equation (4.15) evaluated at the inversion may also be written w' (x,l) = P 1 / 2 Uk 0 h 0 [ (w a ( l ) ) s c a l e d / (w a o (l)) s c a l e d]sin(kx). and results in the same value of the dominant mode as given above. 4.4.3 Particular Example In order to quantify the value of the expressions for the perturbation velocities and the amplitude of the deflection of the inversion, consider parameter values which are again estimated from the model results of C86. For a convective boundary layer depth of 2 km, an L of 1 km, and a nondimensionalized critical wavenumber associated with the roll circulation of 2.221, the corresponding forced wavelength is 5.66 km (0CQ of 1.1105). For a Richardson number of 4.25, the one-layer model results predict a resonant wavelength of 14 km (ocres of 0.448) and a dominant wavelength of 5.74 km ( a d o m of 1.094). Since r|(a) <* w(l) / k the ratio of the amplitude of the forced mode T|0 (a 0 ) and the free modes Ti(a) may be written 1/2 a w « ( D F a d ) 1/2 (4.16) where the subscript denotes the value of alpha for which the amplitude of the velocity and the function F are to be evaluated. Using the notation introduced by (4.8), (4.13) may be written JL "Ho P 1/2 a (wa(l)) scaled ( W « „ ( l ) ) s c a l e d From figure 4.10 this may be further simplified in notation to give Chapter 4: The Stable Upper Layer and ... 166 1/2 'H scaled Cn _ ^ I o J scaled (4.17) For a forced mode with oc0 =1.1105 and an (ri 0) S (, a l e d =0.8506, and a dominant mode with a = 1.094 and an r i ^ j value of 0.8513, the ratio of the warping associated with the dominant mode compared with the forced mode becomes r) / r | 0 = (1.001) (3 1 / 2. Since r) 0 is equivalent to h 0 this gives r| = (1.001)ho • (31 / 2. For an initial forced amplitude of (say) 50 meters, and assuming that the resonant mode distributes energy such that the dominant (and all scales) contains 50% of the energy contained in the forced mode, (4.17) gives a warping of the inversion of 12.5 meters. From (4.15) with a U of 5 m/s, the z-component of the perturbation velocity attains a maximum value of approximately 7 cm/s, along the inversion. It would also be interesting to estimate the vertical perturbation velocity at a height of 2 km above the inversion, in order to compare with observed values of Kuettner et al. (1987) upon which the form of the profiles used in the numerical simulations of C86 were based. In particular, observations from the 12'th of June, 1984 indicate maximum vertical velocities of approximately 3 m/s, while numerical simulations of C86 obtained values of 2 m/s at the end of the model run. Since the assumption as to the amount of energy which is contained in each of the modes is arbitrary it is not meaningful to quantify velocity estimates from the current profiles. However, some insight into the question of which mode is anticipated to be observed at various heights within the troposphere may be inferred from figure 4.9 which illustrates the relative drop in energy with height, for both the resonant and the dominant mode presented in figures 4.6 and 4.7. Figure 4.9 suggests that it is the resonant mode which propagates a larger fraction of it's energy to higher regions of the domain. Thus it is anticipated that the resonant mode will dominant the wave Chapter 4: The Stable Upper Layer and ... 167 field in the upper troposphere, owing to the 'trapping' of shorter wavelengths at lower levels. This restriction on the vertical propagation of shorter wavelength modes is exaggerated by the form of the mean wind profile used in the one-layer model, as compared to numerical simulations. An investigation into the vertical form of the wave field resulting from the more complicated three-layer model (section 4.5) would provide insight into the expected length scale observed at this height. Figure 4.9 further suggests, that if the dominant mode and resonant mode contain equal amounts of energy at the inversion, then the dominant mode will be the observed mode in the lower (approximately) 2 kilometers of the domain. This concept is perhaps more readily applicable to the three-layer model and will be re-examined in section 4.4.2. (Note: The fact that oc0 ~ a d o m is coincidental and is a result of choosing a boundary layer depth of 2 km. If for example, the boundary layer depth is assumed to be 1 km deep, the forced wavelength would correspond to 2.83 kilometers, while the dominant mode and resonant mode would have wavelength values of 5.74 and 14 kilometers respectively. Then a similar calculation as outlined above results in a warping of the inversion associated with the dominant mode of 28.25 meters and a w ' m a x (1) of approximately 16 cm/s. Similarly, for an initial h 0 of 25 meters, the warping is estimated at 14.13 meters, and a w' m a x (1) of 8 cm/s.) Chapter 4: The Stable Upper Layer and . 168 N 5 3 h 2 h ~ i — i 1 — i — i — i — i — i — r J i i i i 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 KE/KEo Figure 4.11: Kinetic energy of the resonant and dominant modes. Plot of KE/KEo where K E 0 is the amount of energy at the inversion. The resonant mode (dashed line: longer wavelength) may be seen to propagate more of its energy to higher elevations, while the dominant mode (solid line: shorter wavelength) is trapped at lower levels. Chapter 4: The Stable Upper Layer and ... 169 4.5 A Three-Layer Model of the Stable Upper Layer. The results of the simple one-layer model presented in section 4.4 are encouraging when compared to the observed wavelength values of Clark et al. (1986). However, using a three-layer model, it is possible to more accurately represent the velocity profile and observed N 2 profiles illustrated in figures 4.2 and 4.3. Due to the added analytical complexity, analysis of the three-layer model (figure 4.10) will focus on determining the resonant modes of the system. If the value of the wavelength of one resonant mode is found to differ only slightly from the results of the one layer model, it will be assumed that the results of section 4.4 adequately represent the dynamics of the stable upper layer. Figure 4.10: Mean wind and N 2 profiles used in the three-layer model. Chapter 4: The Stable Upper Layer and ... 170 As illustrated in figure 4.10, layer 1, ze[0,h,] is characterized by a linearly increasing mean wind U(z) = U( l + z / L , ) , and a constant square of the buoyancy frequency N 2 . Layer 2, z e [ h , , h 2 ] is characterized by a constant mean wind U = U( l + h , / L , ) , and a linearly decreasing stability with height, of the form N 2 (z ) = N 2 ( l - (z - h,) / L 2 ) . Finally, layer 3, z e [ h 2 , ° ° ) is characterized by a linearly increasing mean wind U = U ( l + h, / L l ) ( l + ( z - h 2 ) / L 3 ) , and a constant N 2 = c N 0 , where c> 1. Note that the velocity profile of C86 assumes a constant mean wind for all values of z greater than hi. However, the resulting system of equations will then necessarily require an additional boundary condition (at the top of the domain) in order to fix the resonant modes. Numerically, C86 cap the domain with a 'Rayleigh absorber'. In order to determine the resonant mode of the analytical system, the current study utilizes an increasing mean wind in the uppermost region of the domain, i.e. in layer 3. This may be justified from observations in the lower portion of the stratosphere, which shows an increase in the mean wind with height associated with the stratospheric jets (Holton, 1992, figure 12.3.) Layer 1 If we assume a solution for the vertical component of the perturbation velocity which is stationary with respect to the ground and of the form w' (x, z) = w(z) sin kx , then the general form of the solution for the amplitude w(z) may be written w,(Z) = A Z 1 / 2 K i ( 1 ( a 1 Z ) + B Z 1 / 2 L i ^ ( a 1 Z ) (4.16) Chapter 4: The Stable Upper Layer and ... 171 where K i f l (a ,Z) is defined in section 4.4 and L i ( i (a ,Z) is given by (see Abramowitz and Stegun (1964) and Morgan (1947)) 2 sinh(im) or in integral form as Lj (x) = f cosh(u£)exp(xcos(0))d6- f sin(jo.t)exp(-xcosh(t))dt. sinh(UTi) J 0 J 0 Also, A and B are arbitrary coefficients to be determined, Z = (l + z / L , ) with Z e [1,(1 + 11, / L , ) ] , a, = (k 2 L 2 , ) 1 / 2 , [i = (Ri, - 1 / 4 ) 1 / 2 and the Richardson number for layer 1 is N 2 L 2 given by R i , = . Layer 2 The general form of the solution for w(z) for the case where the mean velocity is constant and the square of the buoyancy frequency decreases linearly with height is given by w2($) = CAi($) + DBi($) (4.17) where A i and B i are Airy functions (see Abramowitz and Stegun (1964)), C and D are coefficients to be determined, and the following notation has been introduced ^ = (Ri 1 / 3 Z - r|), Z = (1 - ( z - h , ) / L 2 ) , TI = c c 2 R r 2 / 3 , a 2 = ( k 2 L 2 2 ) 1 / 2 , and R i 2 = t t 2 / 1 N ? L 2 / t n 2 • Note that Z U (1 + h, / L , ) decrease in value from 1 at z = h, to (1 - (h 2 - h,) / L 2 ) for z = h 2 with (h 2 - h,) < L 2 . The Chapter 4: The Stable Upper Layer and ... 172 corresponding range of £ is given by ( R i 1 / 3 - r | ) at z = h, and ( R i 1 / 3 [ l - ( h 2 - h , ) / L 2 ] - r | ) a t z = h 2 . Layer 3 The upper layer of the three layer model is a layer with linearly increasing mean wind and a constant square of the buoyancy frequency, analogous to layer 1 but with increased stability. Explicitly, U = U( l + h 1 /L,)(1 + ( z - h 2 ) / L 3 ) and N 2 = c N 2 for some o l . The bounded solution for w 3(z) is of the form where F is an arbitrary coefficient, % = (l + ( z - h 2 ) / L 3 ) and thus %e[ l , °° ) , a 3 = ( k 2 L 3 ) 4.5.1 Interface and Boundary Conditions Determination of the arbitrary coefficients A , B, C, D, and F from (4.16) (4.17) and (4.18), requires implementation of the lower boundary condition at z = 0 and as well as the requirement of continuity of the vertical velocity and the pressure across the interfaces at z = h, and z = h 2 . The velocity and pressure conditions at the interface (Drazin and Reid, 1982) are derived in appendix F, and for z = h, are given by w 3(%) = F K i (a 3%) (4.18) y = (Ri, -1 / 4) 1 / 2 and the Richardson number for layer 3 is given by Ri j = ( U d + . ^ / L , ) ) 2 ' w,(h,) = w 2 (h 1 ) - L , H , U dw t(z) + * . ( z ) L i = - L I H I U -dw 2(z) dz dz z = h Chapter 4: The Stable Upper Layer and ... 173 while at z = h 2 we have w 2 (h 2 ) = w 3 (h 2 ) dw 2(z) dw 3(z) z=h, d z z=h 2 dz where H , = (1 + hj / hl). For topography at the lower boundary given by (4.6), the resulting sets of equations may be written in matrix form as Ay = f where y = [ A , B , C , D , F ] T , f = [-k oUh o ,0,0,0,0] T , and the coefficient matrix A is a 5x5 matrix with nonzero elements given in appendix F. 4.5.2 Particular Example Determination of the resonant modes of the three-layer model requires identifying the values of the wavenumber 'k' for which the coefficient matrix A is singular. This was done using values of the parameters estimated from the profiles of section 4.2 and are given in table 4.1. Parameter Layer 1 Layer 2 Layer 3 hi (1km) - - -h 2 (7km) - - -z [0,1] [1,7] [7, infinity) Z [1,2] [1 decreasing to 0.1] — $ — (3.71-3.23k2) to (0.371-3.23k2) — — — [1, infinity) U (m/s) 5(1 + z) 10 10(l + 2(z-7)) N 2 (s-2) 1.15e"4 1.15c-4 ( l - ( z - l ) / 6 . 6 7 ) 3.45e"4 L(km) 1 6.67 1/2 Ri 4.25 51.1 0.8625 Functional Form K a n d L A i and B i K Table 4.1: Parameter values used to calculate the resonant modes. Three-layer model. Chapter 4: The Stable Upper Layer and ... 174 Again using M A P L E , the value of the det(A) was evaluated over a range of wavenumbers, and the values of k for which the determinant attained a null value were noted. For the wavenumber range [0.3, 9.0] numerous zeros were identified. The five zeros corresponding to the largest five wavelengths are associated with values of k between (0.4494,0.4495), (0.9,1.0), (1.3,1.4), (1.8,1.9) and (2.3,2.4). (Note that the resolution of all but the first zero is a result of computational effort associated with narrowing in on the root of the determinant, and not analytic limitations). Note that the roots of the determinant correspond to kj = 0.4494 (and a wavelength of 14 km) and it's approximate harmonics, which are modified slightly towards larger values of k. Some insight into the anticipated behavior may be obtained from results of Wurtele et al. (1987) for a one-layer model (analogous to the model presented in section 4.4) which utilizes values of the Richardson number such that two resonant modes are observed. Results of Wurtele et al. (1987) indicate that the height at which a particular wave is trapped is related to its wavelength with the larger scales propagating to greater heights. Thus it is anticipated that for the current three-layer model it will be the resonant mode closest (in wavenumber space) to the forced mode which will contain the maximum amount of energy. However, since this mode is also anticipated to be trapped at lower heights within the domain, the larger scale resonant mode is expected to be observed as the primary scale in the upper region. 4.6 Summary and Discussion of Results Results of the one-layer model agree well with those of C86 and K87, differing slightly in the wavelength value in the upper tropopause (i.e. model results of 14 km versus 11 km obtained by C86 for a domain which extends 15 km in the vertical, and 13 km for a domain depth of 25 km), Chapter 4: The Stable Upper Layer and ... 175 while the anticipated dominant scale above the inversion show better agreement (5.7 km versus 7.5 km (numerical, C86) and 6 km (observed, K87). The three-layer model also predicts an observed scale in the upper region of the domain of approximately 14 km. It again should be noted that the similarity of the value of the forced wavenumber and the dominant mode of the one-layer model is coincidental, and is a result of the choice of depth for the boundary layer (2 km) as well as the assumption that the forced mode is associated with the critical value of the wavenumber (as predicted by linear theory for the roll circulation). If instead, the boundary layer was 1 km deep, the dominant scale would be approximately twice the forced scale. When compared to the profiles used by C86, the largest source of discrepancy is associated with the assumed value of the stability parameter N 2 . Although care was taken to ensure reasonable choices of all parameters, it is important to obtain an 'accurate' estimate of the value of the Richardson number, as the wavenumber corresponding to the resonant mode of (at least) the one-layer model is sensitive to \i, which is proportional to R i 1 / 2 . Chapter 5 Summary and Conclusions 5.1 Summary Recall that the objective of the current study (section 2.5.4) was to investigate possible mechanisms by which large aspect ratio (i.e. greater than 4) rolls may develop. The particular mechanism of interest is the hypothesized link between boundary layer and tropospheric dynamics (section 2.5.4) and the analytical modeling approach has been to investigate a two-layer model of the troposphere (section 2.5). The convectively driven PBL The effects of incorporating internal gravity wave influences on the nonlinear dynamics of the convectively driven PBL via a modified Rayleigh-Benard convection model (section 3.1) has been investigated using both a weakly nonlinear perturbation analysis (sections 3.3.3 and 3.3.4), and a fully nonlinear truncated analysis (section 3.3.5). The zeroth order equations (associated with classical Rayleigh-Benard convection) have allowed an investigation into two-dimensional nonlinear processes as a possible mechanism of large aspect ratio roll formation (section 3.8.2). In agreement with results of Getling (1983) and Chang and Shirer (1985) upscale energy transfer was observed only for small values of the Prandtl number (i.e. 0.1) for G / G c = 1.004 and G / G c = 1.52. Thus a two-dimensional mechanism is not able to account for laboratory observations of an increase in scales at slightly supercritical values of the Grashof number. 176 Chapter 5: Summary and Conclusions 177 Linear stability analysis (sections 3.6.2 and 3.6.3) and nonlinear numerical results (sections 3.8.2 and 3.8.3) from the first order and truncated equations indicate that the boundary layer will respond with an induced circulation at the imposed wavenumber in association with the corrugated rigid upper lid. Therefore as a minimum, warping of the inversion will result in multiple-scaling including both classical R-B scales as well as scales of motion associated with the imposed wavenumber at the upper surface and its harmonics. Results from the fully nonlinear truncated equations differ from those of the weakly nonlinear perturbation analysis primarily in the predicted modes associated with the imposed wavenumber 'b'. The perturbation analysis indicates a single scale corresponding to 'b' results from the warped upper surface. However, results from the truncated analysis indicate that the forced scale associated with 'b' and its harmonics influence boundary layer dynamics. Specifically, the harmonic of the imposed wavenumber 'b' closest in wavenumber space to the natural development of the flow (as predicted by the zeroth order equations) contains the maximum amount of energy. Explicitly, the dominant wavenumber (a) is estimated to lie within a band of wavenumbers given by ( 2 / 3 < a / a p <4/3) where ap is the preferred wavenumber associated with classical R-B convection (section 3.8.5). A numerical test case, based on atmospherically observed values of the Grashof number with G / G c =22.8 and an eddy Prandtl number of 1, illustrates the ability of the model to produce realistic bulk mean features of the well-mixed convectively driven PBL; i.e. a surface thermal boundary layer and a well-mixed, isothermal layer occupying ~ 70% of the domain. In order to determine the effects of the warped upper surface on roll dynamics, numerical simulations were conducted using both the fully nonlinear zeroth order equations (representing an Chapter 5: Summary and Conclusions 178 isolated atmospheric boundary layer) and the fully nonlinear truncated equations (representing the coupled boundary layer-gravity wave system). Results of the truncated equations suggest that the preferred wavenumber associated with classical R-B convection for G / G c = 22.8 lies in the range [1.5,2.4]. This range of wavenumbers was unexpected, as results from the zeroth order equations suggest a preferred wavenumber of ~3. Also, the truncated system of equations resulted in dominant rolls at wavenumbers (1.6) which are less than the critical wavenumber (2.221) indicating an increase in scale of -25% as a result of the warping at the upper surface. However, forcing at an imposed wavenumber of 1 resulted in a dominant first harmonic associated with a wavenumber of 2. This result suggests that dominant scales which are greater than twice the scale observed at the onset of convection are not realizable by this mechanism. These results agree with those from numerical tests conducted based on values of the parameters associated with laminar convection (section 3.8.4) at marginally supercritical values of the Grashof number, as noted above. The stable upper layer Results from the simple one-layer linear model of the stable upper layer (section 4.4) show good agreement with numerical results of Clark et al. (1986). Velocity scaling based on kinetic energy considerations allowed for the identification of the dominant mode (section 4.4.2) which has the maximum potential for warping of the inversion. Analysis of the three-layer model (section 4.5) identified a number of resonant modes. It is hypothesized that the resonant mode nearest the forced mode will contain the principal amount of energy. Chapter 5: Summary and Conclusions 179 Although the study of the stable upper layer focused on one particular set of atmospheric conditions, the analytical techniques of chapter 4 may be used to investigate a wide range of plausible atmospheric conditions. 5.2 Discussion Results of the current study suggest that an overlying gravity wave field is not able to account for the formation H R V within the convectively driven planetary boundary layer at a wavelength which is more than 50% larger than that predicted by classical R-B convection. Instead, multiple-scaling is predicted where observed scales include 'classical' ones as well as those associated with the imposed wavenumber and its harmonics. Energy contained by the imposed mode has been shown to be 'small' when compared to the energy of the harmonic closest (in wavenumber space) to the preferred scale associated with the classical R-B problem, thus suggesting the existence of only 'weak', larger scales of motion. The freedom of interaction between all scales of motion accounts for the discrepancy between results from the current study and those of Mourad and Brown (1990) (section 2.4). In their study, Mourad and Brown assumed that the harmonics were only forced by the primary wave and results suggest the possibility of a dominant large-scale mode. However, this assumption of the form of harmonic development has been shown to be restrictive and inhibits free flow development. Results from the current study suggest that if the flow is allowed to develop freely and to interact on a wide range of scales, the dominant wavenumber will be limited to a band of wavenumbers centered on of the preferred scale of fluid motion in the absence of external influences, as noted above. Chapter 5: Summary and Conclusions 180 Sang (1991, 1993) developed a linear analytical model of a coupled two-layer system with noted reference to Clark et al. (1984). Discrepancies between the weakly nonlinear perturbation analysis and fully nonlinear truncated analysis suggest that a linear coupling of the P B L and stable upper layer may lead to misleading results, particularly in association with the role of the harmonics and the ability of upper layer dynamics to dominate boundary layer scale selection. Combining results from each layer of the two-layer model of the atmosphere, good agreement is found between the results of the current study and the interpretation of multiple scaling associated with large wavelength cloud streets (as observed by LeMone and Meitin (1984)) based on numerical results presented by Balaji et al., (1993): i.e. the observed large wavelength (-25 km) cloud formations are a result of a combination of important gravity wave and H R V scales, rather than a direct result of H R V with horizontal scales which are comparable to observed cloud band widths. It is important to note that the intention of the present study was to focus on theoretical developments associated with the physical model. Since a number of numerical models exist (e.g. Clark et al. 1986, Balaji et al., 1993) the tendency towards a numerically dominated study has been avoided. In conclusion, although the influence of gravity waves on boundary layer dynamics (as modeled) may result in H R V with aspect ratios which are larger than -3 , this mechanism is not able to account for rolls with aspect ratios greater than -5 . The current model is believed to represent the maximum influence that the gravity waves can exert on the boundary layer, and model results are robust: i.e. wavenumber selection is limited to a band of wavenumbers centered Chapter 5: Summary and Conclusions 181 about the preferred wavenumber associated with classical R-B convection which is independent of the imposed wavenumber and the amplitude of the warping. 5.3 Extensions of Current Work The current study has investigated both atmospheric boundary layer and gravity wave dynamics and has lead to an identification of the effect of an overlying wave field on H R V scale selection; i.e. that of multiple scaling. Discrepancies between the two-dimensional analysis associated with the zeroth order equations and experimental observations of classical R-B convection suggest the extension of the current analytical techniques to three-spatial dimensions in order to determine the importance of dimensional restrictions on flow dynamics for near critical values of the Grashof number. A three-dimensional analysis would also allow for the study of horizontal planform selection. Results from the atmospheric test case suggest the need for further investigation into the trend of wavenumber selection at high values of the Grashof number for both the classical and modified R-B convection models. With respect to atmospherically observed HRV, one weakness of the modified R-B convection model is the assumption of no cross-roll shear. As a first step, the current model was designed to investigate thermally driven H R V circulation and the results obtained allow for a simple estimate of flow behavior. Inclusion of shear would allow for the transfer of energy between the mean flow and the roll circulation and is a logical next step. In addition to the model results of Clark et al. (1984) and observations of Kuettner et al. (1986), an atmospheric test case based on observations of LeMone and Meitin (1984) and Chapter 5: Summary and Conclusions 182 numerical results of Balaji et al. (1993) is also of interest. Balaji et al. (1993) provides sufficient information to allow for an analysis similar to that of chapter 4, which was based on atmospheric profiles of Clark et al. (1984). However, the numerical test cases of Balaji et al. (1993) were conducted using a mean velocity profile which turns with height and scale selection within the boundary layer was interpreted as the Rayleigh mode predicted by a linear theory of convection in a neutral, sheared layer (e.g. Kuo, 1963). Thus, the inclusion of shear into the current model would be required in order to predict the forced mode associated with H R V prior to warping of the inversion (i.e. the forced scale of the stable upper layer). This particular set of numerical results is also of interest as it contains Fourier spectra of the model output for the vertical perturbation velocity within both the P B L and overlying stable layer. With respect to the stable upper layer, the analysis of the three-layer model is limited to determining the wavenumbers associated with the resonant modes. The vertical form of the resonant modes and in particular the longer versus a smaller wavelength modes are of interest in order to investigate the vertical propagation of energy. Also, a nonlinear analysis of the stable upper layer would provide insight into the transfer of energy between tropospheric gravity wave resonant modes and other scales. The partitioning of energy based on kinetic energy assumptions was presented as an alternative (but not a substitute) to a complete nonlinear analysis. Although, a full dynamic coupling of the two-layer model has not been done, it is felt that the rigid upper lid represents the maximum influence that the stable upper layer will exert on boundary layer dynamics. Thus a relaxation of the restrictions associated with the upper bounding surface (i.e. continuity of vertical velocity and pressure) is expected to result in a weaker forcing than has been considered under the current model framework. References Abramowitz M . , and L A . Stegun (1964): Handbook of Mathematical Functions with Formulas, Graphs and Mathematical Tables. U .S. Department of Commerce, National Bureau of Standards, Applied Mathematics Series #55 Angell, J.K., D.H. Pack, and C R . Dickson (1968): A lagrangian study of helical circulations in the planetary boundary layer. J. Atmos. Sci., 25, 707-717. Asai, T. (1966): Cloud bands over the Japan Sea off the Hokuriku District during a cold air outburst. Pap. Meteor. Geophys. (Tokyo), 16, 179-194. Asai, T. (1970a): Three-dimensional features of thermal convection in a plane Couette flow. J. Meteorol. Soc. Japan, 48, 18-29. Asai, T. (1970b): Stability of a plane parallel flow with variable vertical shear and unstable stratification. J. Meteorol. Soc. Japan, 48, 129-139. Asai, T., and I. Nakasuji (1973): On the stability of Ekman boundary layer flow with thermally unstable stratification. J. Meteorol. Soc. Japan, 51, 29-42. Atlas, D., S-H. Chou, and W.P. Byerly (1983): The influence of coastal shape on winter mesoscale air-sea interaction. Mow. Wea. Rev., I l l , 245-252. Atlas, D., B. Walter, S-H. Chou, and P.J. Sheu (1986): The struture of the unstable marine boundary layer viewed by lidar and aircraft observations. J. Atmos. Sci., 43, 1301-1318. Balaji, V . , J.L. Redelsperger, and G.P. Klaassen (1993): Mechanisms for the mesoscale organization of tropical cloud clusters in G A T E phase III. Part I: Shallow cloud bands. J. Atmos. Sci., 50, 3571-3589. Benard, H . (1900): Les tourbillons cellulaires dans une nappe liquid. Rev. Gen. Scienes Pure Appl, 11, 1261-1271, 1309-1328. Boussinesq, J. (1903): Theorie Analytique de Chaleur, Vol 2. Gauthier-Villars, pgs. 154-176. Brown, R.A. (1970): A secondary flow model for the planetary boundary layer. / . Atmos. Sci., 27, 742-757. Brown, R.A. (1972): On the inflection point instability of a stratified Ekman boundary layer. J. Atmos. Sci., 29, 850-859. 183 References 184 Brown, R.A. (1980): Longitudinal instabilities and secondary flows in the planetary boundary layer: A review. Reviews of Geophysics and Space Physics, 18, 683-697. Briimmer, B. (1985): Structure, dynamics and energetics of boundary layer rolls from KONTUR aircraft observations. Beitr. Phys. Atmos., 58, 237-254. Briimmer, B. , B. Rump, and G. Kruspe (1992): A Cold Air Outbreak Near Spitsbergen in Springtime - Boundary-Layer Modification and Cloud Development. Boundary-Layer Meteorol., 61, 13-46. Busse, F.H. (1967): On the stability of two-dimensional convection in a layer heated from below. J. Math. Physics, 46, 140-150. Carruthers, D.J., and C-H Moeng (1987): Waves in the overlying inversion of the convective boundary layer. J. Atmos. Sci., 44, 1801-1808. Chang, H-R, and H.N. Shirer (1984): Transition in shallow convection: An explanation for lateral cell expansion. J. Atmos. Sci., 41, 2334-2346. Chen, M . M . , and J.A. Whitehead (1968): Evolution of two-dimensional periodic Rayleigh convection cells of arbitrary wave-numbers. J. Fluid Mech., 31, 1-15. Chen, Y - Y (1992): Boundary conditions and linear analysis of finite-cell Rayleigh-Benard convection. J.Fluid Mech., 241, 549-585. Chlond, A . (1987): A numerical study of horizontal roll vorticies in neutral and unstable atmospheric boundary layers. Beitr. Phys. Atmosph., 60, 144-169. Chlond, A . (1988): Numerical and analytical studies of diabatic heating effect upon flatness of boundary layer rolls. Beitr. Phys. Atmosph., 61, 312-329. Chlond, A . (1992): Three dimensional simulation of cloud street development during a cold air outbreak. Boundary-Layer Meteorol., 58, 161-200. Christian, T.W., and R . M . Wakimoto (1989): The relationship between radar reflectivities and clouds associated with horizontal roll convection on 8 August 1982. Mon. Wea. Rev., Ill, 1530-1544. Clark, T.L., T. Hauf, and J.P. Kuetner (1986): Convectively forced internal gravity waves: Results from two-dimensional numerical experiments. Quart. J. Roy. Met. Soc, 112, 899-925. Drazin, P.G., and W.H. Reid (1981): Hydrodynamic stability. Cambridge University Press. Ekman, V.W. (1905): On the influlence of the earth's rotation on ocean currents. Arkiv. Math. Astron. Fysilc., 2, 1-54. References 185 Etling, D., and R.A. Brown (1993): Roll vorticies in the planetary boundary layer: A review. Boundary-Layer Meteorol., 65, 215-248. Faller, A.J . (1965): Large eddies in the atmospheric boundary layer and their possible role in the formation of cloud rows. J. Atmos. Sci., 22, 176-184. Faller, A.J . , and R.E. Kaylor (1966): A numerical study of the instability of the laminar Ekman boundary layer. J. Atmos. Sci., 23, 466-480. Fujimura, K. , and J. Mizushima (1988): Truncation effects in a fourier series expansion of nonlinear waves. Phys. Fluids, 31, 2398-2400. Getling, A . V . (1983): Evolution of two-dimensional disturbances in the Rayleigh-Benard problem and their preferred wavenumbers. J. Fluid Mech., 130, 165-186. Haack, T., and H.N. Shirer (1992): Mixed convective cynamic roll vortices and their effects on initial wind and temperature profiles. J. Atmos. Sci., 49, 1181-1201. Haines, D.A. (1982): Horizontal roll vorticies and crown fires. J. Atmos. Sci., 21, 751-762. Hauf, T., and T.L. Clark (1989): Three-dimensional numerical experiments on convectively forced internal gravity waves. Quart. J. Roy. Met. Soc, 115, 309-333. Hein, P.F., and R.A. Brown (1988): Observations of longitudinal roll vorticies during artic cold air outbreaks over open water. Boundary-Layer Meteorol., 45, 177-199. Hildebrand, P.H. (1980): Multiple Dopplar radar observations of P B L structure. Proc. 19th Conf. Radar Meteorology, Miami, Amer. Meteor. Soc, 656-661. Hi l l , G.E. (1968): On the orientation of cloud bands. Tellus, 20, 132-137. Holton, J.R. (1992): An Introduction to Dynamic Meteorology. Academic Press Inc. Holroyd, E.W. (1971): Lake-effect cloud bands as seen from weather satellites. / . Atmos. Sci., 28, 1165-1170. Hurle, D.T.J, E. Jakeman, and E.R. Pike (1967): On the solution of the Benard problem with boundaries of finite conductivity. J. Fluid Mech., 30, 469-475. Jaeckisch, H . (1968): Waveflow above convection streets, l l ' t h OSTIV Congress. Jaeckisch, H . (1972): Synoptic conditions of wave formation above convection streets. 13'th OSTIV Congress. References 186 Jeffreys, H . (1926): The stability of a layer of fluid heated below. Phil. Mag. Ser. 2, 2, 833-844. Kelly, R.D. (1982): A single doppler radar study of horizontal-roll convection in a lake-effect snow storm. / . Atmos. Sci., 39, 1521-1531. Kelly, R.D. (1984): Horizontal roll and boundary layer interrelationships observed over Lake Michigan. J. Atmos. Sci., 41, 1816-1826. Kelly, R.E., and D. Pal (1978): Thermal convection with spatially periodic boundary conditions: resonant wavelength excitation. J. Fluid Mech., 86, 433-456. Konrad, T.G. (1968): The alignment of clear air convective cell. Proc. Intl. Conf. Cloud Physics, Toronto, Amer. Meteor. Soc, 539-543. Koschmieder, E .L. (1969): On the wavelength of convective motions. / . Fluid Mech., 35, 527-530. Krishnamurti, R. (1970a): On the transition to turbulent convection. Part 1. The transition from two- to three- dimensional flow. J. Fluid Mech., 42, 295-307. Krishnamurti, R. (1970b): On the transition to turbulent convection. Part 2. The transition to time-dependent flow. J. Fluid Mech., 42, 309-320. Krishnamurti, R. (1975a): On cellular cloud patterns. Part 1: Mathematical model. J. Atmos. Sci., 32, 1353-1363. Krishnamurti, R. (1975b): On cellular cloud patterns. Part 2: Laboratory model. J. Atmos. Sci., 32, 1364-1372. Krishnamurti, R. (1975c): On cellular cloud patterns. Part 3: Applicability of the mathematical and laboratory models. J. Atmos. Sci., 32, 1373-1383. Kristovich, D.A.R. (1993): Mean circulations of boundary-layer rolls in lake-effect snow storms. Boundary-Layer Meteorol., 63, 293-315. Kropfli, R.A., and N . M . Kohn (1978): Persistent horizontal rolls in the urban mixed layer as revealed by dual-doppler radar. J. Appl. Meteorol., 17, 669-676. Kuettner, J.P. (1959): The band structure of the atmosphere. Tellus, 11, 267-294. Kuettner, J.P. (1971): Cloud bands in the earth's atmosphere: Observations and theory. Tellus, 23, 404-425. References 187 Kuettner, J.P., P.A. Hilderband, and T.L. Clark (1987): Convection waves: Observations of gravity wave systems over convectively active boundary layers. Quart. J. Roy. Met. Soc, 113, 445-467. Kuo, H.L. (1963): Perturbations of plane Couette flow in stratified fluid and origin of cloud streets. Phys. Fluids, 6, 195-211. Langmuir, I. (1938): Surface motion of water induced by wind. Science, 87, 119-123. Laufersweiler, M.J. , and H.N. Shirer (1989): A simple dynamical model of a stratocumulus-topped boundary layer. / . Atmos. Sci., 46, 1133-1153. LeMone, M . A . (1973): The structure and dynamics of horizontal roll vorticies in the planetary boundary layer. J. Atmos. Sci., 30, 1077-1091. LeMone, M . A . , and W.T. Pennell (1976): The relationship of trade wind cumulus distribution to subcloud layer fluxes and structure. Mon. Wea. Rev., 104, 524-539. LeMone, M.A. , and R.J. Meitin (1984): Three examples of fair-weather mesoscale boundary-layer convection in the tropics. Mon. Wea. Rev., 112, 1985-1997. Lipps, F.B., and R.C.J. Somerville (1971): Dynamics of variable wavelength in finite-amplitude Benard convection. Phys. Fluids, 14, 759-765. Lipps, F.B. (1976): Numerical simulation of three-dimensional Benard convection in air. J. Fluid Mech., 75, 113-148. Mahrt, L . (1986): On the shallow motion approximations. J. Atmos. Sci., 43, 1036-1044. Malkus, W.V.R. (1954): The heat transport and spectrum of thermal turbulence. Proc. Roy. Soc. A., 225, 196-212. Malkus, J.S., and H . Reihl (964): Cloud Structure and Distributions over the Tropical Pacific Ocean. University of California Press, 229 pgs. Markson, R. (1975): Atmospheric electrical detection of organized convection. Science, 188, 1171-1177. Mason, P.J., and R.L Sykes (1982): A two-dimensional numerical study of horizontal roll vorticies in an inversion capped boundary layer. Quart. J. Roy. Met. Soc, 108, 801-823. Mason, P.J. (1985): A numerical study of cloud streets in the planetary boundary layer. Boundary-Layer Meteorol, 32, 281-304. References 188 Melfi, S.H., J.D. Spinhirne, S-H Chou, and S.P. Palm (1985): Lidar observations of vertically organized convection in the planetary boundary layer over the ocean. J. Climate Appl. Meteor., 24, 806-821. Miura, Y . (1986): Aspect rastios of longitudinal rolls and convection cells observed during cold air outbreaks. J. Atmos. Sci., 43, 26-39. Morgan, S.P. (1947): Tables of Bessel functions of imaginary order and imaginary arguement. California Institute of Technology, Pasadena. Mourad, P.D., and R.A. Brown (1990): Multiscale large eddy states in weakly stratified planetary boundary layers. J. Atmos. Sci., 47, 414-438. Nield, D.A. (1968): The Rayleigh-Jefferys problem with boundary slab of finite conductivity. J. Fluid Mech., 32, 393-398. O'Brien, J.J. (1970): A note on the vertical structure of the eddy exchange coefficient in the planetary boundary layer. J. Atmos. Sci., 27, 1213-1215. Ostrach, S. (1964): In Theory of Laminar Flows (ed. F.K. Moore), ch 14. Oxford University Press. Pellew, A. , and Southwell, R.V. (1940): On maintained convective motion in a fluid heated from below. Proc. Roy. Soc, Ser A., 176, 312-343. Pitts, D.E., J.T. Lee, J. Fein, Y . Sasaki, K. Wagner, and R. Johnson (1977): Mesoscale cloud features observed from Skylab. Skylab Explores the Earth, NASA-SP-380, pgs. 479-501. Plank, V . G . (1966): Wind conditions in situations of patternform and non-patternform cumulus convection. Tellus, 18, 1-12. Press, W.H. , S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery (1992): Numerical Recipies in Fortran. The Art of Scientific Computing. Second Edition. Cambridge University Press. Priestly, C.H.B. (1962): Width-height ratio of large convection cells. Tellus, 14, 123-124. Proctor, M.R.E. (1981): Planform selection by finite-amplitude thermal convection between poorly conducting slabs. J. Fluid Mech., 113, 469-485. Puhakka, T. and P. Saarikivi (1986): Doppler radar observations of horizontal roll vorticies in Finland. Geophysica, 22, 101-118. Rayleigh, Lord (1880): On the stability, or instability, of certain fluid motions. Proc. London Math. Soc, 11, 57-70. References 189 Rayleigh, Lord (1916): On convection currents in a horizontal layer of fluid when the higher temperature is on the under side. Phil. Mag. Ser. 6, 32, 529-546. Reinking, R.F., R J . Doviak, and R.O. Gilmer (1981): Clear-air roll vorticies and turbulent motions as detected with an airborne gust probe and dual-doppler radar. J. Appl. Meteorol., 20, 678-685. Rothermel, J., and E . M . Agee (1986): A numerical study of atmospheric convective scaling. J. Atmos. Sci., 43, 1185-1197. Sang, J.G. (1991): On formation of convective roll vorticies by internal gravity waves: A theoretical study. Meteor. Atmos. Phys., 46, 15-28. Sang, J.G. (1993): On the dynamics of convection waves. Quart. J. Roy. Met. Soc, 119, 715-732. Sasaki, Y . (1970): Influences of thermal boundary layer on atmospheric cellular convection. J. Meteorol. Soc. Japan, 48, 492-501. Schluter, A. , D. Lortz, and F. Busse (1965): On the stability of steady finite amplitude convection. J. Fluid Mech., 23, 129-143. Shirer, H.N. , and J.A. Dutton (1979): The branching hierarchy of multiple solutions in a model of moist convection. J. Atmos. Sci., 36, 1705-1721. Stensrud, D.J., and H.N. Shirer (1988): Development of boundary layer rolls from dynamic instabilities. J. Atmos. Sci., 45, 1007-1019. Stull, R.B. (1993): An Introduction to Boundary Layer Meteorology. Kluwer Academic Publishers. Sykes, R.L, W.S. Lewellen, and D.S. Henn (1988): A numerical study of the development of cloud-street spacing. J. Atmos. Sci., 45, 2556-2569. Sykes, R.L, and D.S. Henn (1988): On the numerical computation of two-dimensional convective flow. J. Atmos. Sci, 45, 1961-1964. Walter, B .A. (1980): Wintertime observations of roll clouds over the Bering Sea. Mon. Wea. Rev., 108, 2024-2031. Walter, B .A. , and J.E. Overland (1984): Observations of longitudinal rolls in a near neutral atmosphere. Mon. Wea. Rev., 112, 200-208. Watson, A. , and G. Poots (1971): The effects of sinusoidal protrusions on laminar free convection between vertical walls. J. Fluid Mech., 49, 33-48. References 190 Weston, K.J. (1980): An observational study of convective cloud streets. Tellus, 32, 433-438. White, D.B. (1988): The planforms and onset of convection with a temperature-dependent viscosity. J. Fluid Mech., 191, 247-286. Willis, G.E., J.W. Deardorff, and R.C.J. Somerville (1972): Roll-diameter dependence in Rayleigh convection and its effect upon the heat flux. J.Fluid Mech., 54, 351-367. Woodcock, A . H . (1942): Soaring over the open sea. Sci. Mon., 55, 226-232. Wurtele, M.G . , R.D. Sharman, and T.L. Keller (1987): Analysis and simulations of a trophosphere-stratosphere gravity wave model. Part 1. J. Atmos. Sci., 44, 3269-3281. Appendix A Scaling Based on Atmospheric Parameters of Surface Heat Flux and the Convective Velocity Scale w*. Recalling that the primary motivation for the study of the modified R-B convection problem was to parallel phenomenon observed within the convectively driven PBL, consideration is given to the choice of scaling based on atmospherically derived parameters such as the surface heat flux (w'0') s and the convective velocity scale (w*). A . l Surface Heat Flux For the convective boundary layer the surface heat flux may be parameterized using where U is the magnitude of the mean horizontal wind at a height of z above the surface (where z is usually taken to be 10 meters), C H is a bulk transfer coefficient (with typical values ranging from l x l O - 3 to 5xlO - 3 ) , 6 is the mean potential temperature (assumed constant with height), and 9 S is the potential temperature of the surface. (See Stull, 1988.) However, for the problem under consideration, the dimensional surface heat flux may be calculated explicitly and is assumed given by where 5Q denotes the surface heat flux, T is the temperature field, K is the thermal diffusivity (either eddy or molecular thermal diffusivity, depending on the problem under consideration), fi is 191 ( w e ' ) s = - c H u ( e - e s ) (A.1) 8Q = - K V T n (A.2) Appendix A: Scaling Based on Atmospheric ... 192 the unit normal pointing into the fluid, the temperature gradient is evaluated at the lower plate, and the bar denotes spatially averaged. A.2 Velocity Scale The atmospheric convective velocity scale may be written w* = [ ( g / 9)(w' 8 ' ) s Z j ] where z{ is the depth of the boundary layer and (w' G')s is given by (A.l) . (See Stull, 1988.) For a surface heat flux 5Q given by (A.2) the convective velocity scale is assumed of the form 7^(oQ)rI 1/3 (A.3) where T m has replaced the mean potential temperature of the C B L , and the depth of the boundary layer z s is naturally replaced by the mean distance H between the bounding surfaces. A.3 Estimate of the Convective Velocity Scale for Rayleigh-Benard Convection Using the experimental data presented by Willis et al. (1972), assume that the external parameters are given by H = 2.54xl0~ 2m, A T m a x = 3 0 ° C (for all experiments conducted), and a Tmean =300°K (room temperature). Further assume that a = 0.71, v = 1.461xl0"5m2 / s , and K = 2.06xl0" 5m 2 /s where a is the Prandtl number (v / K ) , the kinematic molecular viscosity is denoted by v , and K is the molecular thermal diffusivity. Based on these values of the parameters the velocity scale based on a viscous dissipative time scale of -41 seconds is approximately b x l O ^ m / s , or O(10"4) where the 'O ' denotes 'is of order'. Using (A.2) the surface heat flux based on a AT of 0.5°C, is Appendix A: Scaling Based on Atmospheric ... 193 8Q ~ 4 . 0 6 X 1 0 " 4 m°C / s and the convective velocity scale becomes w* = 7 x l 0 - 3 m / s, or O(10" 3), with a convective time scale of 3.5 seconds. A .4 Estimate of the Convective Velocity Scale for Atmospherically Observed H R V For an eddy viscosity of 100 m2/s (Miura, 1986) and a boundary layer depth of 1 kilometer, the velocity scale is estimated to be 0.1 m/s with a corresponding time scale of 2.8 hours. For a surface heat flux of approximately 300 watts/m2, and a mean potential temperature of 300 K, w* = [(g / 6)(w' 0' ) s Zj ] 1 / 3 over the same boundary layer depth gives a value of w* = 2 m / s, with a corresponding time scale of 8.3 minutes. (Note that the velocity scale based on the (eddy) viscous dissipative time scale is an order of magnitude less than that obtained for w*.) A.5 Nondimensionalization Based on the length scale H and the convective velocity scale (w*) the following nondimensional variables are introduced X = * , Z = ± , 0 = 2 ^ R x = ^ t , . *F = ^ - V . H H ( T 0 - T , ) H ( w * ) H Y The resulting nondimensionalized form of the governing equations for the streamfunction and temperature are (compare with (3.11) and (3.12)) J ( ¥ , V 2 ¥ ) + £ 2 G ^ - i W 4 v F = 0 (A.4) dx dX — - J O F , 0 ) - - V 2 0 = O (A.5) dx a where J is the Jacobian, G is the Grashof number and r> as the ratio [v / [(w*)H]]. Appendix B Equations Governing the Real and Imaginary Components of the Zeroth Order, First Order, and Truncated Equations B . l Zeroth Order Equations Letting 9 a k = 9 R ° + i9 a ° k and \ | / a k = \ | / R ° + h|/ a° k in equations (3.31) and (3.32), leads to the following for the real and imaginary components of the temperature and streamfunction: 1} ^ r - = + 2 a V U - - X X J[ cUn,m , k K-,) , m e R , n +v ( R a-a , , m eI , n ]da ' m=l n=l Z\Q I I 1 00 co so 2) ^ = — Y . 2 X k - 2 a ¥ a R k +^XX J[cU..„.^IvS-..xmel% - v U ^ e ^ J i a ' m=l n=l . 3) - Y l k ^ L = aG9:, k +Ya4kVf,k jY^„[cU,m,k][V(Ra-o,m¥l,n +¥U),m¥aR,„]da' m=l n-1 . 4 ) - Y a . k = _ a G e a . k + Ya.kVl.k ~ ~ X X j Y ^ [Cly,n,m,k Ivfa-a'Xm Va^.n ~ V j a - a ' X m V a ' . n ] d a ' (B.l) where y 2 q = (p 2 + ( q 7 t ) 2 ) , y 4 k = ( y 2 k ) 2 , 1 1 ^ n m k = Jcos(m7rri)sin(n7rri)sin(k7rri)dri, ( j ) ~ m k = Jsin(m7tr|)cos(n7rri)sin(k7cn)dr|, 0 0 and C a , a , n > m , k =[ (a-a ' ) (n7t ) ( t ) ; m , k -a ' (m7t)< m , k ] . (B.2) 194 Appendix B: Real and Imaginary Components . B.2 First Order Equations Similarly, equations (3.35) and (3.36) become r)AR 1 1 ! ) - T ~ = T a , k 0 a , k +2av|/a,k dt cr 71 X X JC l.a',n,m,k[V(Ra-a'),m 0a> + V ' a - a ' X m 6 ^ " + V f a - a ' X m 0 " , , + ¥ |a-a' ),m 0 a\n ] d a ' + 2 h R ( V » 9 " ) where h R ( < k , 9 ° k ) = - b 2 a 2 0 _ [ 5 ( a _ b ) + 5( a + b)] +-(k7t)2[e RO 4_QR 0 (a-b),k D(a+b),k +^X(n7t)(3n,k [(a - b/2)9R0_b),n - (a + b/2)9R 0 + b ), n ] " n=l 1 V 1 V 1 f [ p 2 - RO +C2+ i i f R 0 fa10 Ha' • / > / i J |_ a,a',n,m,k T (a-a'-b),m a,a',n,m,k T (a-a' + b),m _p a'.n m =l n=l _oo L V Y f fc 2 - \ i / 1 0 + C 2 + \ i / 1 0 FJR Oda' J L a,a',n,m,k T(a-a'-b),m a,a',n,m,k T (a-a'+b),m j^a'.n m=l n=l 2) = Y a A . k -2a^/ a > k 9T CT +iy y fc1 , rvf eR0 -v!1 e10 + \ i / R 0 e R 1 - w ! 0 e1! Ida' J a,a ,n,m,k [_ Y (a-a ),m a ,n Y (a-a ),m a ,n Y (a-a ),m a ,n Y (a-a ),m a ,n J ft m=l n=l . + 2 h , « k , e ° k ) where h R « k , 9 ° k ) = - ( Ioc ) 2 [ e » b M l + 9| a° + b ), k] +-X(mt)PBik [(a - b / 2 ) 9 » b ) , n - (a + b/2)9| a° + b ), n ] G n=l Appendix B: Real and Imaginary Components 196 + — - Y Y [\c2~ \ i / R 0 + C 2 + v / R 0 "bR Oda' A Z-lZ-i) L a , a ' , n , m , k Y ( a - a ' - b ) , m a ,a ' ,n,m,k T ( a - a ' + b ) , m _ P a \ n ^ft m = l n = l _ „ , XX J [ ^ a , a ' , n , m , k V ( a - a ' - b ) , m + ^ a , a ' , n , i i i , k V(a-a'+b),ra ^ a ' . n ^ 3 ^ft m = l n = l _o 3 ) - y 2 M L a , k dx a G e i : k + y : , k ¥ R k + i y y i y , rC' , ir v R i „ v 1 0 + X I / J 1 v i / R 0 Ida' ' • ' . J I a ,n L a,a , n , m , k J ^ T ( a - a ) , m Y a ,n Y ( a - a ) , m Y a ,n J ft m = l n = l + ^ XX jya2,n[Cl,,,n,m,kIvra-,),mV^ + ¥ ! a 0 - , ) , m ¥ a R , , „ ] d a ' +2q R « k , K k ) 7 1 m = l n = l -where qR«k,e;,k) = -^2(n*>MeS-». " 6 ^ ] - ( k 7 t ) 2 y a 2 k [ V ( R O _ b ) , k + < + b ) , k ] 2 tr, +X(n7t)Pn,k[a-VR0_b),„ + ~ 2 2 J [ « 3 ¥ ( R ° - , - b ) , m +a> ( R a°- , + b), m]¥^da' n = l ^ft m = l n = l _o= T ^ X X J [«3¥:L -b) , r a + « : ¥ ! a ° - , + b ) , m ] ¥ a R ° „ d a ' m-1 n = l . « k 4) - Y l k ^ = - a G e R , k + y : , k ¥ i : k _ i y y i y , [c1 , T v R l V R 0 - vi/!1 V° Ida' J I a , n |_ a,a , n , m , k JL T ( a - a ) , m Y a ,n Y ( a - a ) , m Y a ,n J ft m = l n = l 1 X X j Y a . , n [ c I , , , n , m , k ] [ V ( R 0 - , ) , m ¥ a R , 1 „ " V:a°-a,,mVI>]da' + 2 q ' « k , 0 ° k ) ft m = l n=r where q I « k , e ° k ) = Gb 2 [ 5 ( a + b) - 8(a - b)] + 2(™)Pn,k \[e J°_bXn - e(Ra°+b),n ] Appendix B: Real and Imaginary Components ... 197 - (kTt ) 2 Ya2k [<_b),k + < + b ) , k ] + X ( ™ ) P , k [«2 V j - b X . + «2 ¥!a+b),„ ] n=l + ^ X I JfcKVb),™ +OC>(Ra°-,+b),m]¥aROndaI '^'t m=l n=l (B.3) m=l n= The coefficients not defined by (B.2) are given by: p\ k = Jricos(n7tri)sin(k7rn)dr| o I I X k , n , m = Jtlsin(n7rn)sin(m7rTi)sin(k7ni)dri X k , „ , m = Jtlcos(n7cri)cos(m7cri)sin(k7cn)dri with C 2 ; a , n > m , k = ( m 7 t ) ( a ' ) ^ ; m , k -(nTt)(a-a'-b)(t)-,m | k 2^+ ''a,a',n,m,k = (m7t)(a' )(t>;m k - (im)(a - a'+b)^ m.k and a2 - • 2b(a - b )y ( V B ) ; „ - 2b j (a - b) - b 2 (n r t ) 2 - 3b z(a - b) a ; - — + 2b(a + b)YU,n +2bJ(a + b ) - b z (n7 t ) z -3b 2 (a + b)J 1 (a')y2 -(a'+bXnrc)2 +<l>;m . k (a -a ' -b ) (n7c) iy2,n+(a')b + (n7c) 2+ib 2 +Xn.m,k (a-a'-b ) (n7u) 2 -b (a ' ) - -b 2 2 + X ; m . k ( n 7 t ) -b(a ' )2 (m7t ) - | (a ' )b 2 - ^ b 3 ( i r w T ) « 3 = <t>;,m,k - - ( a ' )(nur)Y2 n - (a' )(m7t)(n7t)2 +b(m7c)(n7c)2 Appendix B: Real and Imaginary Components . 198 + Y ; m , k ( a - a ' +b ) (n7u) lYa 2 , „ - ( a ' ) b + ( n 7 t ) 2 + i b 2 +Xtm,k(a-a'+b)(n7i:) : ( a ' ) b - - b 2 2 +X;m, k(n7t) b(a' )2 (irax) - | b 2 (a") + ^  b 3 (nm) (B.4) B.3 Complete Equations Finally, equations (3.37) and (3.38) become: D ^ L = - 1 Y a 2 k e a R k + 2 a < k dT CT XX j[CU,„,m,kl¥!a-a),m0R,„ +¥(Ra-a),me:,„]da' + £ 2h R ( ¥ a k , 0 a k ) m=l n=l . rift1 1 2) -r— = —Ya,k9a,k - 2a\|/a,k dx CT J[cU,„,m,k][vR-a),m0aR,„ - vU),m0l,„]da' +e2h1(ij/a,k,ea,k) ^ m=l n=l . -T \ „ , 2 ^ a , k „ , ^ Q I , , , 4 , „ R 3 ) - Y a , k — ^ — = aGe a k + Y a . k V a . k dX + - 2 Z lYa 2 , „ [C l ,a . , „ , m , k ] [¥ ( R a-a ' ) , M V^,n + Y | . - . > V . . , + ^ (V . .k • ^ a . k ) ^ m=l n=l 4) - Y 2 , k ^ = - a G e R , k + Y a 4 k ¥ I , k j Y a , n [ C l y , n , m , k ] [ ¥ ( R a - a ' ) , r „ V R n - V U x r n V l . n J d a ' +£2q I (V|/AK , 0 a k ) ^ m=l n=l where h R ( ¥ a k , 0 a k ) , h ' ty^ .e^) , qR0|/ a,k> f ia, k) and q'(vi/ a i k ,0 a , k) are defined in sectionB.2. Appendix C Simplifying — (V + e V ) ^ arising from the Truncated Equations dx In order to solve for the Fourier coefficients using explicit methods, the term — ( V 2 dx resulting from the mapping of equation (3.15), was simplified by neglecting eV*P. Such a simplification may be justified by considering ^ 2 a2 * a 2 ^ which in wavenumber space is approximated as ( a 2 +(k7t ) 2 ) l \ | / k l where l \ | / k l denotes the magnitude of the first mode of the streamfunction. Now d2xV d2xP ?W V ¥ = 2nbsin(b<;)4—-2cos(bc;)^-T- + Tib 2 C os (b$)^ -and since lcos(b£)l< 1, isin(b^)l < 1, and 0 < r\ < 1 we have that ™T, ^d2x? ^d2V u 2 a ^ w < 2 , W V + b v T E R M (1) TERM(2) T E R M (3) Consider T E R M (1), which may be approximated in wavenumber space as 2 t > l ^ n ~ 2 a b ( k 7 t ) l x t / k l • Similarly T E R M (2) becomes 2-|-^- < 2(k7U) z l \ | / k l and finally, T E R M 2i a y (3) reduces to b 2 — < b 2 (k7t ) l \ | / k l , 199 Appendix C: Simplifying ... 200 From the results presented in section 3.8, we find that |\|/,| » |\)/> 1| (where the subscript ' > 1 ' denotes vertical modes greater than one) assume l\)/ kl corresponds to l y , I . Then the requirement that V 0 2 v F » eV,*?, leads to (a 2 + 7 U 2 ) I \ | / 1 I » £7r,[2ab + 27t + b 2 ] \ | / 1 l . Thus the restriction on the amplitude of the forcing becomes ( a 2 +7C 2 ) 8 « 7c[2ab + 27t + b 2 ] which for 0 < a < 10 and 0 < b < 2 gives 8 « 0.063. Appendix D The Unstable Conductive Base State associated with the First Order Equations The inabil i ty o f the f l ow to support a horizontal temperature gradient may be readily seen by returning to the original f o r m o f the governing equations given by (3.1) through (3.4) . Consider the x and z components o f the momen tum equations d u 3t • + 3u du U h W dx dz aw lit •+ aw aw dx dz + — ^ - v V 2 w + ( p P m ) j _ 3 p P m ^ g = o and the equation o f state g iven by ( p - p m ) _ ( T - T J I n order that a stable conduct ive solut ion exists, i.e. u = w = 0 , requires aP dx o - ^ = - ( p - P m ) g (D.l) (D.2) D i f fe rent ia t ing (D.2) w i t h respect to x, apply ing ( D . l ) , and substi tut ing f r o m the equation o f state, leads to the condi t ion dx Conversely, i f a hor izontal temperature gradient exists w i th in the f l o w , the x component o f the m o m e n t u m equations states that the resultant hor izontal pressure gradient must be balanced by a nonzero hor izontal component o f veloci ty (i.e. u ^ 0 ) . 201 Appendix D: The Unstable Conductive ... 202 This result may be compared by analogy to flow between vertical plates of uniform temperature (e.g. Ostrach, 1964) by considering the vertical temperature gradient along line A and line B as shown in figure D. 1. B T o Figure D . l : Horizontal temperature gradient associated with the warped upper plate Since the temperatures along the upper and lower surfaces are assumed constant, and Ah | A < Ah| B where h is the distance separating the two plates, the resultant vertical temperature gradients are such that 9 T dz 3T dz (D.3) Consider points 1 and 2 located at a height d above the lower surface, on lines A and B respectively, as illustrated in figure D.2. 1 • * 2 Figure D.2: Horizontal temperature at a height'd' above the lower plate. Appendix D: The Unstable Conductive ... 203 As a result of (D.3), T p o i n t l < T p o i n t 2 and therefore there exists a horizontal temperature gradient along the plane joining points 1 and 2. D . l Comparison with Experimental Conditions A natural question arises as to the applicability of the above results since experimental limitations include geometric imperfections which according to the above analysis, would result in an induced circulation regardless of the magnitude of the horizontal temperature gradient. However, the above analysis also assumes that the plates are perfectly conducting. Thus, the degree of geometric imperfections that may be supported is related to the degree of imperfection in the conductivity of the bounding surfaces. A conductive solution (which requires horizontal isothermal planes) may then be realized, as illustrated in figure D.3. Figure D.3: Experimental environment. To estimate the necessary limitation on the conductivity of the plate such that a stable conductive solution may exist, consider the steady-state temperature equation (3.3) given by K V 2 T = 0 which (for small amplitude of forcing) along the surface of the plate may be written Appendix D: The Unstable Conductive ... 204 9 2 T a 2T n where the subscripts 'p' and 'a' denote plate and air respectively. Rewriting, this becomes K p A T b 2 = - K a ^ where (1/b) is the imposed wavelength. Thus ( a 4 ) Example Consider an experimental scenario with water between aluminum plates. From the CRC Handbook of Chemistry and Physics 59 t h ed., 1978, the thermal conductivity of water is 6.09 (mW) / (cm°K) while the thermal conductivity of aluminum is 2.37 W / (cm°K). Thus for an imposed wavenumber b = 1 / cm (D.4) gives e< -^-S—=0.51 mm. Therefore, deflections of 0.51 mm over 1 cm would still support a conductive solution. Appendix E Comparison of Analytical and Numerical Linear Results This appendix presents the results of the comparison between analytically obtained solutions to the first order and truncated equations with the numerical results. E . l First Order Equations For comparison of numerically obtained values of the coefficients of temperature and the streamfunction with the analytical results, both the conductive and convective zeroth order solutions were considered. Conductive zeroth order base state For the assumptions outlined in section 3.6.2, with G = l , a = 1, b = 1 and Cla initial conditions, the analytical and numerical results (from tests F l and F120) are summarized in tables E . l and E.2. Recall that in the conductive regime, the only nonzero component of the forcing (f,) occurs in association with the imposed wavenumber 'b', as given by (3.49). Similarly, results for G = 250, a = 1, b = 1 and Cla initial conditions, (from tests F2, F200 and F210) are summarized in tables E.3 and E.4. Both cases show excellent agreement between numerical and analytical results. It should be noted that the existence of a nonzero value of the forcing at 'b' occurs within one time step. In physical space, this represents the simultaneous, and instantaneous, formation of a uniform roll field over the entire domain. Also, in comparison to the rate at which an unforced flow attains it's 205 Appendix E: Comparison of Analytical... 206 'steady-state' value, the forced circulation stabilizes rapidly. (Compare for example, the classical nonlinear R-B convection model for a Grashof number of 1.004 times the critical value, which has not yet reached steady-state, after a nondimensional time of 200.) y components Analytical f. Numerical f. 1 -4 -4.000 2 0 0.000 3 0 0.000 4 4 4.000 Table E . l : Analytical and numerical components of the forcing vector f,. Parameter values: G = l , o" = 1, b = l and Cla. Variable Analytical Numerical T = 1 x final 1 T =2 1 final -0.3748 -0.3748 -0.3748 e b 0 0.25e"5 0.92e"10 ¥ b 0 -0.10e"5 -0.25e 1 0 . -0.0370 -0.0370 -0.0370 Table E.2: Comparison of analytical and numerical results for a = b . Parameter values: G = 1, a = 1, b = 1 and Cla. Appendix E: Comparison of Analytical. 207 y components Analytical f, Numerical f. 1 -4 -4.000 2 0 0.000 3 0 0.000 4 1000 1000.0 Table E.3: Analytical and numerical components of the forcing vector f,. Parameter values: G = 250, c = 1, b = 1 and Cla. Variable Analytical T =1 1 f ina l 1 Numerical T =2 1 final ^ T -4 1 final ^ e bR -3.153 -3.091 -3.152 -3.153 e b 0 0.44e"3 0.74e"5 0.21e8 v R 0 -0.15e"2 -0.25e"4 -0.7 l e 8 -15.136 -14.926 -15.132 -15.135 Table E.4: Comparison of analytical and numerical results for a = b . Parameter values: G = 250, a = 1, b = 1 and Cla. Convective zeroth order base state For the assumptions outlined in section 3.6.2, with G = 330, a = 1, b = 1 and the amplitude of the zeroth order temperature and streamfunction coefficients, as determined from test N30, as A R = A 1 = 0.04525 and B R = ( -B 1 ) = 0.1511. The analytical and numerical results for the first Appendix E: Comparison of Analytical... 2 0 8 order variables, using Cf l l initial conditions are summarized in tables E.5 and E.6 through E.8. As with the conductive test, analytical and numerical values of the solution vectors agree. f. Analytical f3 I f, Numerical I I 1 -4 0.4160 0.4658 -4.000 0.4160 0.4658 2 0 0.4160 0.4658 0.000 0.4160 0.4658 3 0 28.305 16.670 0.000 28.305 16.670 4 1320 -28.305 -16.670 1320.000 -28.305 -16.670 Table E.5: Analytical and numerical components of the forcing vectors f(, f , f4. Parameter values: G = 330, a = 1, b = 1 and CIII. Variable Analytical '''final — 2-Numerical ^ final = 4. tfina. = 6 -eR -4.986 -4.975 -4.986 -4.986 el 0 0.000 0.000 0.000 0 0.000 0.000 0.000 -25.100 -25.054 -25.100 -25.100 Table E.6: Comparison of analytical and numerical results for a = b . Parameter values: G = 330, a = 1, b = 1 and CIII. Appendix E: Comparison of Analytical. 209 Variable Analytical Numerical ^ f ina l — 2- ' ' ' f i n a l = 4- t f i n a l = 6 -0 ( R a „ - b ) 0.4014 0.3961 0.4014 0.4014 0 ( a o - b ) 0.4014 0.3961 0.4014 0.4014 K - o ) -1.5035 -1.4830 -1.5032 -1.5035 < 0 - b ) 1.5035 1.4830 1.5032 1.5035 Table E.7: Comparison of analytical and numerical results for (a + b) = a 0 . Parameter values: G = 330, o" = 1, b = 1 and CIII. Variable Analytical Numerical ^ final — 2- ^ final = 4- X final = 6 " Q R ° ( a „ + b ) 0.5086 0.4943 0.5082 0.5086 D ( a „ + b ) 0.5086 0.4943 0.5082 0.5086 < „ + b ) -1.4682 -1.4271 -1.4670 -1.4682 < 0 + b ) 1.4682 1.4271 1.4670 1.4682 Table E.8: Comparison of analytical and numerical results for (a - b) = a 0 . Parameter values: G = 330, a = 1, b = 1 and CIII. Appendix E: Comparison of Analytical... E.2 Truncated Equations For the special case of m = 3, the coefficient matrix B of (3.48) becomes 210 B = Yo -eA7 1 lb 0 0 0 eGb/4 0 0 E [ A ; | O + A - | O ] 0 0 Y 2 b - e A " | 0 0 Gb eGb/4 0 - e A r Y 2 b - e A , 0 -eA+| 12b 3b 0 - e G b / 4 yi o o 0 0 0 -yt eB~ 0 -2ab 0 0 e K - B « U -Y: 0 0 0 .0 -4ab 0 0 -6ab 0 0 eBr G2b - e G b / 4 0 eGb / 4 G3b 0 eB; 0 12b ' Y 2 b 0 e B : i I 2b eB: 13b Y^b with y 2 = [ e R , e R , 0 R b , 0 R b , \ ) / o , \ | / b , \ ) / 2 b , \ | / 3 b ] T and f, =[0,-4b2e,0,0,0,4Gbe,0,0]T. Conductive regime Using a horizontal wavenumber cutoff of 3, analytical results (section 3.6.3) based on (3.46) and (3.47), and a comparison with numerically obtained values of the Fourier coefficients associated with 'b' and it's harmonics, are presented. For the assumptions as outlined in section 3.6.3, with G = l , a = 1, N = l , b = 1, e = 0.01 and T f i n a l = 4 , the analytical and numerical results from test C l using Cla initial conditions are summarized in table E.9. Similarly, results for G = 250, rj = l , N = l , b = l , £ = 0.01, and T f i n a l =8, from test C2 using Cla initial conditions are summarized in table E.10. Appendix E: Comparison of Analytical Variable Analytical Numerical e0R -1.52c-4 -1.52e"4 eR -3.75e"3 -3.75e-3 oR -5.30e"5 -5.30e"5 oR -5.37e'7 -5.37e"7 ¥o ¥,' 0 -3.70e"4 0 -3.70e"4 ¥ 2 -5.5 le"6 -5.5 le"6 vl -5.62e"8 -5.63e"8 Table E.9: Comparison of analytical and numerical results for a = 0, b, 2b, 3b. Parameter values: G = 1, a = 1, b = 1, £ = 0.01, T f = 4 and Cla. Variable Analytical Numerical e0R of oR ¥0 v! ¥ 2 ¥ 3 -1.29e"3 -3.18e"2 -4.2 le"3 -2.49e"4 0 -1.52c"1 -1.31e"2 -6.54e"4 -1.29e" -3.18e": -4.21e": -2.49e"' 0 -1.52e" -1.31e": -6.55e"' Table E.10: Comparison of analytical and numerical results for a = 0, b, 2b, 3b. Parameter values: G = 250, a = 1, b = 1, £ = 0.01, T f = 8 and Cla. Appendix E: Comparison of Analytical... 212 Convective regime For values of the parameters given by G = 330 , a = 1, N = 1, b = 1, 6 = 0.01 and T f i n a l = 12 , the analytical and numerical results from test C30 using Class III initial conditions are summarized in table E. 11. Variable Analytical Numerical e 0R or e 2R e R ¥ o V i ¥ 2 V3 -3.19e"j -7.86e"2 -4.63c"1 -8.28e2 0 -3.6061 -1.60 -2.46e_1 -2.3 le" -5.70e -1.17e" -1.96e": 0 -2.78e" -4.02e" -5.82e": - 2 Table E . l 1: Comparison of analytical and numerical results for a = 0, b, 2b, 3b. Parameter values: G = 330, a = 1, b = 1, e = 0.01, T f = 12 and CIII. Appendix F Derivation of the Interface Conditions, and the Elements of the Coefficient Matrix from the Three-Layer Model Test Case F . l Derivation of the Interface Conditions Recall the x-component of the steady-state momentum equation, as well as the equation of continuity given by (4.2) as - d u ' , dU 1 dp' A U — + w ' — + — ^ r - = 0 (F.l) ox dz p 0 dx ^ + ^  = 0 (F.2) dx dz For w' (x,z) = w(z)sin(kx), the equation of continuity (F.2) suggests a horizontal component of the perturbation velocity of the form u'(x,z) = u(z)cos(kx), and from (F.l) a perturbation pressure dependence given by p' (x,z) = p(z)cos(kx). Explicitly, (F.2) suggests that the vertical dependence of the perturbation velocity components are related by ^ = ku ( z ) . (F . 3 ) dz Substitution of the pressure, and vertical velocities into (F.l) and simplifying gives 3 T J 1 -Uku(z) + w(z) — = — kp(z). (F.4) Po Using (F.3) to replace u(z) in (F.4) leaves =dw(z) „ dU 1 - U — — + w ( z ) — = —kp(z) . (F.5) dz dz p 0 213 Appendix F: Derivation of the Interface 214 The condition of the continuity of pressure across an interface thus becomes - dw^z) dU, - U i — z — + w,(z)-dz dz z = h - dw 2(z) „ - U 2 — + w 2 (z) dz d U 2 dz (F.6) z = h where the subscript ' 1' denotes values of the flow variables evaluated at the interface (at z = h ) from below, while the subscript '2' denotes values of the flow variables evaluated at the interface from above. The second interface condition requiring continuity of vertical velocity is [ w l ( z ) ] z = h = [ w 2 ( z ) ] ^ h . (F.7) Specific form of the interface conditions for the three-layer model The interface conditions given by (F.6) and (F.7) are now applied to the interfaces located at z = hj and z = h 2 of the three-layer model of the stable layer as outlined in section 4.5 and depicted in figure 4.9. Recall that the general form of the solutions for each of the three layers is given by w 1 (Z) = A Z , / 2 K i t l ( a 1 Z ) + B Z 1 / 2 L i ^ ( a 1 Z ) (F.8) w2($) = CAi($) + DBi£ ) w 3(x) = F K i Y ( a 3 x ) (F.9) (F.10) where the parameters are defined following (4.16), (4.17) and (4.18), and A, B, C, D, and F are arbitrary coefficients to be determined. • Lower interface at z = hj. From (F.8) and (F.9), application of (F.7) gives Appendix F: Derivation of the Interface .. 215 A H 1 1 / 2 K i ( i ( a , H 1 ) + B H 1 1 / 2 L i t l ( a 1 H 1 ) = CAi(l) + Bi(l) /2 (F . l l ) where since Z = (1 + z / Lj) the value of Z at z = h, is given by Hj = (1 + z / L , ) . Similarly £ = 1 when z = h, . When applied to the lower interface (F.6) becomes -U(l + h, / L , ) — A 1 i / T , , , 3K ( (a,Z) - H r 1 / a K ^ ( a 1 H 1 ) + H , l / a ' ^ z ' Z = H , -U(l + h, / L , ) — B | L , ^ H r , / a L i , ( a I H 1 ) + H 1 I " * 9Z Z = H , + ^ [ A H 1 1 / 2 K i / a 1 H 1 ) + B H / / 2 L ^ (a.H,)] 1/3 ua+h^L,)^- ,3Ai($) + D 3Bi(§) 5=i 5=i (F.12) • Upper interface at z = h 2 . From (F.9) and (F.10) application of (F.7) gives C A i ( H 2 ) + DBi (H 2 ) = F K i y ( a 3 ) (F.13) where H 2 is the value of ^ | z = h , and we note that % = 1 at the upper interface where z = h 2 . Finally, application of (F.6) gives = F[(l / 2)K i y (a 3 ) - K i Y (a 3)] / L 3 (F. 14) Ri 1/3 ,3Ai(0 + D 3Bi(0 5 = H 2 5 = H 2 Appendix F: Derivation of the Interface ... 216 Lower boundary condition Application of the lower boundary condition given by (4.6) gives A K ^ ( a 1 ) + B L i | l ( a 1 ) = - U h 0 k (F.15) F.2 Elements of the Coefficient Matrix The resulting sets of equations may be written a„A + a 1 2 B = g (F.16) a 2 1 A + a 2 2 B + a 2 3 C + a 2 4 D = 0 (F.17) a 3 1 A + a 3 2 B + a 3 3 C + a 3 4 D = 0 (F.18) a 4 3 C + a 4 4 D + a 4 5 F = 0 (F.19) a 5 3 C + a 5 4 D + a 5 5 F = 0 (F.20) where g = - U k h 0 , or in matrix form as Ay = f (F.21) where y = [ A , B , C , D , F ] T , f = [-kUh 0,0,0,0,0] T, and the coefficient matrix A is a 5x5 matrix with nonzero elements given by (F.12) through (F.15) as a„ .= K i ( 4 ( a I ) a 1 2 = L ^ ( a , ) a 2 1 = H j / 2 K i ( i ( a 1 H 1 ) a 2 2 = H ^ L ^ H , ) a 2 3 =-Ai (^ ( l ) ) a 2 4 =-Bi(^( l ) ) a 3 1 = - H , [(1 / 2 ) H : , / 2 K i ( l ( a , H , ) + H | / 2 K ' ^ (o^rl,)]/ L , +[H 1 1 / 2 K i ^ (a,!!,)]/ L , a 3 2 = - H 1 [ ( l / 2 ) H - 1 / 2 L i / a 1 H , ) + H ; / 2 L , i J l ( a 1 H 1 ) ] / L 1 + [ H , I / 2 L i H ( a 1 H 1 ) ] / L 1 Appendix F: Derivation of the Interface ... 217 a 3 3 = - H . R i f A f G(D) / L 2 a 3 4 = - H . R i f B i ' £ (1 ) ) / L 2 a 4 3 = Ai(£(H 2 )) & 4 4 = Bi($(H 2)) a 4 5 = - K i y ( a 3 ) a 5 3 = R i 2 / 3 A i ' ( ^ ( H 2 ) ) / L 2 a 5 4 = R i 2 / 3 B i ' ( £ ( H 2 ) ) / L 2 a 5 5 = [ ( l / 2 ) K i Y ( a 3 ) - K ' i y ( a 3 ) ] / L 3 where the prime denotes differentiation with respect to the function's argument (i.e. d I dZ for K and L , and d I dt, for the Airy functions). Also H , = (1 + h, / L j ) , £(1) = [Ri 2 3 - r |], and ^ (H 2 ) = [ R i 2 / 3 ( l - ( h 2 - h , ) / L 2 ) - T i ] . 

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:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0087910/manifest

Comment

Related Items