A C O N V E C T I V E B O U N D A R Y L A Y E R M O D E L W I T H T I M E - A D A P T I V E GRIDS By M A G D A L E N A R U C K E R B. Math, University of Waterloo, 1992 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTERS OF SCIENCE in THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF MATHEMATICS, INSTITUTE OF APPLIED MATHEMATICS We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA January 1995 © M A G D A L E N A R U C K E R , 1995 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, 1 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 The University of British Columbia Vancouver, Canada DE-6 (2/88) Abstract This thesis presents the first implementation of time-adaptive non-uniform grid dis-tributions in higher order turbulence closure modelling of the convective boundary layer. The aim was to introduce a vertical grid that resolves the surface and entrainment layers, but uses coarser resolution elsewhere. Two types of grid generating methods were exam-ined with the Level 2.5 turbulence closure model (Mellor and Yamada, 1982; Yamada, 1983). The solution profiles of mean and turbulent quantities were compared to model results obtained with stationary grid distributions of the same mesh size as well as a reference solution obtained with a grid resolution of 10m. ii Table of Contents Abstract ii List of Tables v List of Figures vi List of Abbreviations and Major Symbols ix Acknowledgement x 1 Introduction 1 1.1 The Convective Boundary Layer 2 1.2 Rationale 5 1.3 Research Objective and Outline 6 2 Model Description 7 2.1 Governing Equations 9 2.2 Numerical Scheme and Boundary Conditions 12 2.3 Model Performance 14 3 Mesh Selection Methods 17 3.1 A Priori Mesh Selection 20 3.2 Adaptive Mesh Selection 21 3.2.1 Theoretical Background 21 3.2.2 Basic Implementation 23 iii 3.2.3 Weight Function 24 4 Results 26 4.1 Invariant Grids 26 4.2 Adaptive Mesh Method 31 4.3 A Priori Mesh Method 39 5 Summary of Results and Conclusions 47 Bibliography 49 iv List of Tables 4.1 Hourly C B L height for log-linear grid distributions 28 4.2 Hourly C B L height for adaptive grids with N = 44 34 4.3 Hourly C B L height for adaptive grids with N — 22 38 4.4 Hourly C B L height for a priori grids, Case I 42 v List of Figures 1.1 Typical vertical structure of (a) mean potential temperature 0 and (b) turbulent sensible heat flux w9 in a C B L . w90 is the surface heat flux; w9i is the entrainment heat flux; Z; is the C B L height.(Source: Cai, 1993) . . 3 2.2 Location of variables on staggered grid 12 2.3 Vertical profiles for (a) mean wind component U and (b) Reynolds stress component uw 15 2.4 Effects of location of first grid point z\ and grid resolution Az on u*. The horizontal line represents the value of obtained with a log-linear resolution 16 4.5 Log-linear grid distributions for Az^ = 10, 55 and 129m corresponding to c = 0.2, 0.0355 and 0.015, respectively 27 4.6 Vertical profiles of mean wind velocity components (a) U and (b) V ; and Reynolds stress components (c) uw and (d) vw for grid distributions shown in figure 4.5 29 4.7 Vertical profiles of (a) mean potential temperature 0; and (b) vertical heat flux w6 for grid distributions shown in figure 4.5 30 4.8 Vertical profiles of T K E for grid distributions shown in figure 4.5 30 4.9 Grid distributions for adaptive method with N — 44, a = 0.0001 and a — 0.0005. Corresponding log-linear grid as well as reference profiles are included 32 vi 4.10 Vertical profiles of (a) mean potential temperature 0; and (b) vertical heat flux w9 for grid distributions shown in figure 4.9 32 4.11 Vertical profiles of mean wind velocity components (a) U and (b) V; and Reynolds stress components (c) uw and (d) vw for grid distributions shown in figure 4.9 33 4.12 Grid evolution for adaptive method with N — 44, (a) a — 0.0001 and (b) a = 0.0005 34 4.13 Grid distributions for adaptive method with N = 22, a = 0.0001 and a = 0.0005. Corresponding log-linear grid as well as reference profiles are included 35 4.14 Vertical profiles of (a) mean potential temperature 0; and (b) vertical heat flux w9 for grid distributions shown in figure 4.13 36 4.15 Vertical profiles of mean wind velocity components (a) U and (b) V; and 'Reynolds stress components (c) uw and (d) vw for grid distributions shown in figure 4.13 37 4.16 Grid evolution for adaptive method with N = 22, (a) a = 0.0001 and (b) a = 0.0005 38 4.17 Grid distributions for a priori method, Case I. fTg is the regridding fre-quency; cubic interpolation is indicated by cu and quadratic interpolation by qu. Corresponding log-linear grid as well as reference profiles are included. 40 4.18 Vertical profiles of mean wind velocity components (a) U and (b) V; and Reynolds stress components (c) uw and (d) vw for grid distributions shown in figure 4.17 41 4.19 Vertical profiles of (a) mean potential temperature 0; and (b) vertical heat flux wd for grid distributions shown in figure 4.17 42 vii 4.20 Grid distributions for a priori method, Case II. fTg is the regridding fre-quency; cubic interpolation is indicated by cu and quadratic interpolation by qu. Corresponding log-linear grid as well as reference profiles are included. 44 4.21 Vertical profiles of mean wind velocity components (a) U and (b) V; and Reynolds stress components (c) uw and (d) vwior grid distributions shown in figure 4.20 45 4.22 Vertical profiles of (a) mean potential temperature 0; and (b) vertical heat flux wB for grid distributions shown in figure 4.20 46 viii List of Abbreviations and Major Symbols A B L Atmospheric boundary layer C B L Convective boundary layer LES Large eddy simulation T K E Turbulent kinetic energy L Monin-Obukhov length N total number of nodes or mesh size q*/2 turbulent kinetic energy u x-component of mean wind velocity ug x-component of geostrophic wind uw x-component of Reynolds stress V y-component of mean wind velocity vg y-component of geostrophic wind vw y-component of Reynolds stress w6 turbulent sensible heat flux WQQ turbulent sensible heat flux at the surface w9i turbulent sensible heat flux at the inversion base friction velocity Zi height of CBL, or height of inversion base roughness length Z\ location of first mean grid point Az vertical grid spacing Az^ vertical grid spacing for linear section of a log-linear grid 0 mean potential temperature K von Karman constant ix Acknowledgement I wish to express my gratitude to a number of people who have contributed advice, support and effort to this research, academically or in other ways. I am very grateful to my research supervisor, Dr. Douw G. Steyn, for his guidance, inspiring enthusiasm and financial support in form of Research Assistantships throughout my program. I would also like to thank Dr. Brian Wetton for taking time to serve on my supervisory committee. His suggestions were greatly appreciated. Many thanks are extended to Dr. Xiaoming Cai who was always available with an attentive ear and offered invaluable advise on numerous occasions. Special thanks to all for the prompt reading of this thesis. During the course of this study, the Department of Mathematics provided personal funding through several Teaching Assistantships. Dr. Keith Ayotte at the National Center for Atmospheric Research, Boulder, Col-orado, provided the numerical code for the Level 2.5 turbulence closure model which formed the basis of my work. I thank him for his helpfullness and generosity of his time and knowledge. Many thanks must also go to Vincent Kujala who provided much assis-tance in numerical and computational matters and who finally made the printing of this thesis possible. Lastly, I would like to extend my deepest gratitude to Allen Fujimoto for his contin-uous support and encouragement. x Chapter 1 Introduction The atmospheric boundary layer (ABL) signifies that portion of the atmosphere which is influenced by frictional and thermal effects of the earth's surface. Almost all biological life occurs in this layer and is therefore directly affected by its processes. A good under-standing of the development and structure of the A B L is essential to many fields such as agricultural and air pollution meteorology, weather forecasting and climate studies. The development and structure of A B L has been subject of an enormous number of observational (Stull, 1988), analytical (Tennekes and Driedonks, 1981) and numeri-cal investigations (e.g., Yamada and Mellor, 1975; Andre et al., 1978; DeardorfF, 1972). Much of the earlier work concentrated on the characteristics of turbulence close to the surface because measurements were readily attainable there (Caughey, 1982). As interest grew in the turbulence dynamics and structure of the upper parts of the A B L , so did the use of numerical models such as higher order turbulence closure models and Large Eddy Simulation (LES) models. In higher order turbulence closure models, the Reynolds stresses which result from applying Reynolds decomposition and ensemble-averaging to the Navier-Stokes equations are solved prognostically. The equations for the Reynolds stresses, however, contain unknown third order moments; equations for the third order moments contain unknown fourth order moments and so on. This dilemma is known as the closure problem (Stull, 1988) and is solved by parameterizing the unknown variables in terms of some combination of known variables. In volume-averaged LES models, the large scale turbulence fields are resolved explicitly on a three-dimensional grid, while the 1 Chapter 1. Introduction 2 subgrid-scale effects are parameterized (Zeman, 1981). Higher order turbulence closure models have enjoyed much popularity in A B L modelling as they provide information on the vertical distribution of turbulent statistics (Garratt, 1992), and are computationally faster than LES models (Zeman, 1981). This study concerns itself with the effects of variable grid resolution in the numerical modelling of a specific state of A B L , called the convective boundary layer (CBL). The following section provides further background on the C B L . Sections 1.2 and 1.3 outline the rationale and objective for this research, respectively. 1.1 The Convective Boundary Layer The structure of the A B L is strongly influenced by features of the dominant mechanism driving the turbulence. Turbulence is created mechanically through wind shear or ther-mally through the upward transfer of heat. On sunny days, when the ground is strongly heated by solar radiation, convective mixing dominates mechanical mixing. Under this condition, the A B L is commingly referred to as the C B L . The typical structure of the C B L can be represented by vertical profiles of mean potential temperature and turbulent sensible heat flux (figure 1.1). Large convective plumes, generated by surface heating, extend vertically over the whole C B L , thus causing intense vertical mixing. As a result, the mean potential temperature and mean wind velocity profiles are nearly independent of height; this part of the C B L is referred to as the mixed layer. Almost all the wind shear and potential temperature gradients are confined to a shallow region close to the ground, called the surface layer. In the entrainment layer, jumps in mean wind velocity and potential temperature are apparent as they adjust from the height independent values in the mixed layer to values characteristic of the free atmosphere. Chapter 1. Introduction 3 (a) Mean / \ T Temperature / J Inversion r / \ Layer (b) Turbulent Sensible Heat Flux ^ Entrainment i i Layer \ $ __ _ L \ Mixed w'6 ; Layer \ Surface Layer r \ \ Figure 1.1: Typical vertical structure of (a) mean potential temperature 0 and (b) tur-bulent sensible heat flux u>9 in a C B L . w0o is the surface heat flux; w6j is the entrainment heat flux; Z; is the C B L height.(Source: Cai, 1993) , Chapter 1. Introduction 4 The mixed layer is often capped by a stable layer or inversion which limits the vertical extent of convection and therefore defines the height of the C B L , Zi. The C B L depth exhibits strong diurnal variation, ranging from a few tens of metres in the early morning to 1 - 2 km by mid afternoon on a sunny day. The growth of C B L is controlled by the surface sensible heat flux w9o, mean subsidence and the stability of the air aloft. The growth of the C B L is associated with the entrainment process, whereby convective elements continuously penetrate into the stable layer aloft, pulling warm, turbulence free air down into the cooler, mixed layer. The entrainment of warm air into the cooler air below results in a negative (downward) entrainment heat flux below the inversion layer. The ratio of entrained to surface heat flux is approximately —0.2 for convective conditions (Driedonks and Tennekes, 1984). The region which is dominated by small scale entrainment effects is called the entrain-ment or interfacial layer and extends approximately from 0.8.Z; to 1.2Z,- (Caughey, 1982 ). Higher order turbulence statistics, however, show entrainment effects deep down into the mixed layer (Driedonks and Tennekes, 1984). This downward extension of entrainment effects on higher order turbulent statistics is still considered as unresolved (Driedonks and Tennekes, 1984). It should be emphasized that the properties and structure of the C B L are influenced by surface forcings as well as by entrainment. A positive surface heat flux as well as a negative entrainment heat flux will increase the mean potential temperature in the mixed layer. A negative surface momentum flux together with a positive entrainment momentum flux will change the mean wind velocity in the mixed layer. Chapter 1. Introduction 5 1.2 Rationale In numerical modelling, generally speaking, the accuracy of solution depends on two constituents: the physical model and the numerical model. The physical model entails the parameterizations and assumptions which define and justify the system of differential equations to be solved. The numerical model involves the discretization of the physical model on a discrete grid. The grid must be sufficiently fine so that the discrete problem represents the continuous problem with adequate accuracy. The size of the grid directly influences the computational cost and thus, it is desirable to choose a mesh which trades off accuracy with cost. For solutions with narrow regions of high spatial activity, an equispaced fixed mesh with high resolution is computationally inefficient, and depending on the complexity of the physical model and duration of simulation, not always feasible. A more computationally efficient method would be to locally refine the grid. Moreover, if regions of high spatial activity are moving in time, then the grid also needs to adapt in time. In C B L modelling, two regions of high spatial activity exist: the surface and entrain-ment layer. While the former region is fixed in time, the position of the latter is time dependent; that is, it moves upward during the daytime as air from the free atmosphere is entrained into the mixed layer. As pointed out in the previous section, the processes at the surface as well as in the entrainment layer affect the development and structure of the C B L . Thus, one may expect it necessary to represent both the surface and entrainment layer with sufficiently fine resolution in order to achieve an accurate solution. The need for sufficiently fine grid resolution is amplified for higher order turbulence closure models for which the closure depends on local gradients. In previous numerical studies of the C B L , fixed log-linear grid distributions have been employed to accurately resolve the gradients at the surface (e.g., Andren, 1990; Chapter 1. Introduction 6 Yamada and Mellor, 1975). The non-stationarity of the entrainment layer, however, has not allowed a similar implementation of local grid refinement without resorting to a uniformly spaced fine mesh. Thus, the effects of grid resolution in entrainment modelling remain a concern, especially for higher order turbulence closure models (Ayotte, Andren; pri vate communi cat ion). 1.3 R e s e a r c h O b j e c t i v e a n d O u t l i n e The objective of this study is to investigate the effects of non-uniform grids in higher order turbulence closure modelling of the C B L . Two different types of grid generating methods are examined that produce fine resolution in the surface and entrainment layers and coarser resolution elsewhere. The Level 2.5 turbulence closure model (Mellor and Yamada, 1982; Yamada, 1983) is used for simplicity. The model is one dimensional and does not include effects of humidity or radiation. The study will proceed as follows. The turbulence closure model used in this study, as well as details of the numerical scheme and boundary conditions are described in Chapter 2. Chapter 3 discusses two grid selection methods. Results of the grid selection methods as well as comparisons to uniform grid distributions are contained in Chapter 4. Finally, a summary of results and conclusions are given in Chapter 5. Chapter 2 Model Description In higher order turbulence closure models, closure is achieved by parameterizing the unknown moments in terms of the lower moments and their gradients. Thus, in a Ith turbulence closure technique, the (J + l)th moments are parameterized in terms of Ith or lower moments or their local gradients. It is assumed that as the closure approximations are of higher orders, the equations for lower order variables become more accurate (Lum-ley and Khajeh-Nouri, 1974). Although in theory, closure can occur at an arbitrarily high order, certain limitations must be considered in practice: • For / > 3, the turbulence closure models are far too complex to be used in large scale applications. For example, the 3 r d order turbulence closure model for a moist atmosphere by Andre et al. (1978) involves 33 differential equations for turbulence variables. • Accurate observational measurements become increasingly difficult with higher or-der moments. Due to the lack of physical observations, closure approximations for higher orders may be crude and physically meaningless (Stull, 1988). • The higher the closure, the more parameters there are to "tune" the model. For these reasons, the second order closure technique has received the most attention. It is the lowest order that provides information on higher order turbulence statistics, but still includes parameterizations based on physics. 7 Chapter 2. Model Description 8 Various levels of approximations and simplifications exist even for second-order tur-bulence closure modelling. Mellor and Yamada (1974) (1982) developed a hierarchy of second-order turbulence closure models based on the isotropic dissipation hypothesis by Kolmogoroff (1941) and the energy distribution hypothesis by Rotta (1951). In this hi-erarchy, the full second order turbulence closure, which is termed Level 4 and consists of ten prognostic equations1 is systematically simplified using anisotropy scaling. Level 3 neglects the time rate of change, advection and diffusion terms for all but the isotropic components of turbulent moments. Only the turbulent kinetic energy (TKE) and po-tential temperature variance are computed prognostically, while anisotropic components are calculated algebraically. For Level 2, time rate of change, advection and diffusion terms are dropped even for the isotropic components, which reduces the entire turbu-lence model to an algebraic set of equations. Mellor and Yamada (1982) introduced a Level 2.5 model, which retains a prognostic equation for T K E , but neglects the time rate of change and the vertical diffusion of the potential temperature variance. The Level 2.5 model has enjoyed much popularity as it represents a balanced compro-mise between accuracy and numerical efficiency (Helfand and Labraga, 1988). Its appli-cations include environmental studies (Yamada, 1977), air pollution modelling (Andren, 1990) and general circulation models (Miyakoda and Sirutis, 1977). However, numer-ous deficiencies have been reported with this closure. For example, Zeman and Lumley (1976) point out that the simple gradient transport assumption for third order fluxes is inadequate for buoyancy driven mixed layers. Helfand and Labraga (1988) show that the Level 2.5 model produces non-realizable turbulent fields2 in case of rapidly growing turbulence. x This assumes a dry atmosphere. For moist atmosphere, the number of prognostic equations is higher. 2 The realizability condition requires that the variances must always be non-negative and all fluxes must vanish at the same time that T K E vanishes. Chapter 2. Model Description 9 For this study, the Level 2.5 closure model was chosen to investigate the effects of grid resolution as it represents a well recognized and relatively simple closure model. Flaws associated with this model are irrelevant, as comparisons with other models or observational measurements are not undertaken here. It is assumed that the findings of this study with respect to grid resolution may be applicable to other higher order turbulence closure models. 2.1 Governing Equations The governing equations for the mean variables in a dry, horizontally homogeneous at-mosphere with no subsidence are dU ,.T, . duw s i - f i y - v.) - -gp (2.D where upper case symbols represent mean quantities and lower case fluctuating quantities. Thus, U and V are the mean horizontal winds, 0 is the mean potential temperature, uw and vw are the Reynolds stresses, and w6 is the vertical heat flux. / represents the Coriolis parameter. Horizontal pressure gradients have been replaced by the geostrophic wind components (Ug, Vg). For simplicity, the x-axis is aligned with Ug. Radiative effects are neglected. For the Level 2.5 model, a prognostic equation for T K E is added dq2 dU dV q3 d = -2uw— - 2vw— + 2g/3w9 - 2 f - + — dt dz ^ dz, ^ — ' ^jV^ dz 5, dq2 (2.4) Chapter 2. Model Description 10 where q2 = u2 + v2 + w2 represents twice the T K E , g is the gravitational acceleration, (3 — 1/0 is the thermal expansion coefficient, and A x and X1 are turbulent length scales. The right hand side of (2.4) consists of shear production (I), buoyant production or loss (II) and frictional dissipation (III). The diffusion term (IV) represents redistribution by transport and pressure forces. The diagnostic equations for the velocity variances, momentum fluxes, heat fluxes and potential temperature variance can be algebraically manipulated to obtain expressions for the Reynolds stresses and vertical heat flux in equations (2.1) to (2.3) uw = -KM-Q^ (2.5) dV vw = - K M — (2.6) ^ B = - K H ^ - . (2.7) dz The eddy diffusivities KM and KH are defined by KM = lqSM (2.8) KH = lqSH. (2.9) where q is the square root of twice the T K E and / is the master turbulent length scale which will be defined later. This study uses the version of Level 2.5 model after Yamada (1983), whereby the sta-bility functions SM and SH are substituted from the Level 2 model in order to circumvent the problem of non-realizable turbulence fields. SM and SH are algebraically expressed in terms of the flux Richardson number Rj S h = (2-10) 1-R f Chapter 2. Model Description 11 b M " A2 5 l 7 i - [Si(7i + 72) - 3 ]^^ b H where 7 1 = I - (2A1/B1) (2.12) 7 2 = (B2/B1) + (QAi/Bi). (2.13) The flux Richardson number Rf is derived from the gradient Richardson number Ri through the relation Ri = ~Rf (2-14) where = gpae/dz . 15) The critical flux Richardson number occurs when SM = SH = 0 in equations (2.10) and (2.11) and is obtained byi?/c = 71/(71 +72)- Thus, Rjc = 0.191 defines the Richard-son number above which turbulence and mixing cease to exist. Following Yamada (1983), Rf is kept below 0.16 in order to provide background values for SM and SH-The length scales Ai, A2, h, h and Ai are taken to be proportional to a master turbulent length scale / so that (/1; Ai, 12, A2, \i) = (A\, B^ A2, B2, E{)1. A commonly used length scale formula is that of Blackadar (1962) which interpolates between the two limits / ~ KZ as z —+ 0 and I ~ /0 as z —> 00 where K is the von Karman constant and the asymptotic value lo is given by lo = 0 . 2 ^ . (2.17) Jo <ldz Chapter 2. Model Description 12 •t + 1 q2, uw, vw, wO, KM, KH / / / / / / / / / / Figure 2.2: Location of variables on staggered grid. The constants of proportionality are determined emperically from neutral data. For this study, values from the Mellor and Yamada (1982) paper are assumed: (Ai,Bu A2, B2, Ei) = (0.92,16.6,0.74,10.1,0.01). The constant C\ is taken as 0.08. 2 . 2 N u m e r i c a l S c h e m e a n d B o u n d a r y C o n d i t i o n s A staggered grid is used with mean quantities calculated at main levels and turbulent quantities at intermediate levels (figure 2.2). The first mean grid level is positioned at height z\. The surface Reynolds stresses and vertical heat flux are defined at height ZQ, where ZQ is the roughness length. Central differencing is used to compute spatial derivatives while t ime integration occurs wi th an implici t forward Euler scheme. In the equations for mean quantities, the eddy diffusivities are taken from the previous time step. This allows for a tridiagonal matrix which can be efficiently inverted. A semi-implicit centered scheme with weight Chapter 2. Model Description 13 0.75 on the future time step is used for the diffusion and dissipation terms in the T K E equation. Because a semi-implicit scheme for the dissipation can cause negative T K E , a check is performed to ensure that T K E is positive at all times 3. Initial vertical profiles for mean variables are taken from a L E S simulation of a C B L and some modifications were made: the inversion base was shifted from a height of approximately 1200m to 460m to prevent the inversion base from reaching the upper boundary during a simulation. The lapse rate was set to a value of 0 .003° K / k m . T K E is initialized using the algebraic Level 2 turbulence closure model (Yamada, 1983), in which the production of T K E balances the dissipation. The model is driven by a constant surface heat flux of W9Q = 0 . 0 8 ° K m / s and barotropic geostrophic forcing of (Ug, Vg) = (10,0) m/s. Due to the presence of explicit terms in the discretization, the scheme is not unconditionally stable. However, a conservative time step of 5 seconds was found adequate for the resolutions used. Length of simulation is 5 hours for all cases in this study. Upper boundary conditions for the model are: U = U„ V = Vg, ^ = 0, q2 = 0. (2.18) oz At the lower boundary, the surface heat flux is specified as input. The friction velocity u* is obtained iteratively using mean wind velocities at z\ and Monin-Obukhov similarity theory through the relation (Paulson, 1970) K (2.19) where, for convective conditions, * „ = -21o (1±£) - In (l±f!) + 2 t a „ - , - \ 3i .e. 9— > e where £ is a small positive number. Here, e = 1 0 - 5 . Chapter 2. Model Description 14 a n d x - I 1 - 15 L. \\U\\ is t h e n o r m o f m e a n w i n d v e l o c i t y a n d L is t h e M o n i n - O b u k h o v l e n g t h d e f i n e d b y L = —u^/(k^rwd0). T h e s u r f a c e R e y n o l d s s t resses a re d e t e r m i n e d f r o m t h e f r i c t i o n v e l o c i t y u * u s i n g UWQ = — ul cos a i , vw~o = — ul s i n a i (2 .20) w h e r e Q i is t h e a n g l e b e t w e e n t h e su r f a c e a n d g e o s t r o p h i c w i n d d i r e c t i o n s , d e t e r m i n e d at t h e f i r s t g r i d p o i n t z\. B o u n d a r y c o n d i t i o n for T K E at t h e l o w e r b o u n d a r y is s p e c i f i e d t h r o u g h t h e r e l a t i o n fo r c o n v e c t i v e c o n d i t i o n s g i v e n b y A n d r e n (1990 ) f ul ^ = (1 - 3z/L)2'3 + 2 (6 .2 - 0 . 5 2 Z 4 / I ) 2 / 3 (2 .21) 2.3 Model Performance I n t h e s u r f a c e l a y e r , i t is a s s u m e d t h a t t h e R e y n o l d s s t resses a r e a p p r o x i m a t e l y c o n s t a n t ; t h a t i s , t h e y v a r y b y less t h a n 1 0 % . T h u s z\ is o f t e n a s s i g n e d a v a l u e w h i c h l i e s i n t h e s u r f a c e l a y e r , b u t i s o t h e r w i s e a r b i t r a r y 4 . T h e l o c a t i o n o f z\ a n d t h e g r i d r e s o l u t i o n Az, h o w e v e r , h a v e b e e n f o u n d t o af fect t h e R e y n o l d s s t resses as w e l l as m e a n w i n d v e l o c i t y c o m p o n e n t s . F i g u r e 2.3 c o m p a r e s t h e v e r t i c a l p r o f i l e s o f t h e m e a n h o r i z o n t a l w i n d c o m p o n e n t U a n d t h e R e y n o l d s s t ress c o m p o n e n t uw fo r u n i f o r m g r i d d i s t r i b u t i o n s o f 1 0 m a n d 1 0 0 m . F o r z^ = 5 m , b o t h v e r t i c a l p r o f i l e s o f U a n d uw w i t h c o a r s e r r e s o l u t i o n di f fer c o n s i d e r a b l y f r o m t h e p r o f i l e s w i t h 1 0 m r e s o l u t i o n . S i m i l a r d i s c r e p a n c i e s c a u s e d t h r o u g h t h e u s e o f 4 eg . Z l = Chapter 2. Model Description 15 4 6 8 10 12 U (m/s) -0.4 -0.3 -0.2 -0.1 0.0 0.1 uw (m2/s2) Figure 2.3: Vertical profiles for (a) mean wind component U and (b) Reynolds stress component uw. an arbitrary mesh near the lower boundary were also pointed-out by Taylor and Delage (1971). To alleviate this problem, they recommend the use of a constant flux wall layer (for neutral case) or fine grid spacing near, the boundary as found with a log-linear grid. Indeed, Figure 2.3 shows almost identical profiles for 10m resolution and a log-linear resolution where the grid spacing increases to a constant value of approximately 60m. Figure 2.3 also shows that the differences between profiles with different grid resolution greatly decrease if z\ is moved further away from the boundary. This observation is further illustrated in Figure 2.4, which shows the friction velocity for various Z\ and Az. Inaccurate values of u* are obtained with coarse grid resolution near the surface combined with a small z\. This discrepency becomes larger with increasing Az and decreasing Zi, but virtually disappears for large values of z\. Chapter 2. Model Description 16 © X o A X © x 0 • A • dz = 10 A dz - 20 X dz = 50 O dz = 100 6 g g r— 1 1 1 — " 1 0 10 20 30 40 50 2, (m) Figure 2.4: Effects of location of first grid point z\ and grid resolution Az on u„. The horizontal line represents the value of obtained with a log-linear resolution. This study is mainly concerned with the effects of grid resolution above the surface layer. For comparison purposes, it is desirable to keep z\ constant. Thus, the effects of grid resolution and location of z\ in the surface layer, as demonstrated above, are removed by comparing meshes that are logarithmically spaced near the surface, but approach different linear grid spacing away from the surface. Chapter 3 Mesh Selection Methods An optimal mesh can be theoretically defined as a mesh which yields a minimum mesh size N for a given error tolerance TOL. The exact solution to this optimization problem, however, is not feasible (Ascher et al., 1988) and so instead we look for a good mesh which is not much larger than the optimal grid size N for a given TOL. In practice, a slightly different, but related formulation of the mesh selection problem is of more use: for a given mesh size N, find a mesh which minimizes the error. It has been observed that the choice of such a mesh is often not unique in that, for a given problem and TOL, a wide range of acceptable meshes of size N may exist, even when a uniform mesh of the same size yields poor results (Ascher et al., 1988). Mesh selection methods can be grouped into two categories: a priori and adaptive. If the spatial and temporal behaviour of the solution is known, a mesh may be chosen a priori which appropriately conforms to the expected behaviour of the solution. Difficul-ties in finding an appropriate transformation may arise, however, if the solution is badly behaved (Ascher et al., 1988). An example of a priori mesh selection in numerical C B L modelling is the inclusion of a time independent, log-linear grid distribution in order to resolve strong surface gradients (e.g., Andre et al., 1978; Yamada and Mellor, 1975). In adaptive methods, the physics of the problem at hand control and direct the grid distri-bution so that the physical solution is represented with sufficient accuracy. In general, some solution quantity related to the error is equidistributed1 in order to minimize the 1 Equidistribution may be considered as the unifying feature for adaptive methods in general (Thomp-son, 1985). 17 Chapter 3. Mesh Selection Methods 18 error. The concept of error equidistribution originated from boundary-value problems in ODE's (Thompson, 1985) and is described by Ascher et al. (1988). Adaptive methods may be divided into two further subclasses: static-regridding and dynamic-regridding methods. In static-regridding methods, the grid is adapted at dis-crete time steps and interpolation is required to transfer numerical quantities from the old to the new grid. In dynamic-regridding methods, also known as moving mesh methods (Zegeling, 1991), the discretization of PDE's and the grid selection procedure are intrin-sically coupled. The grid points move continuously in the course of the solution in order to properly resolve developing solution gradients or higher derivatives, and interpolation is not required. These two adaptive approaches have their own advantages and drawbacks. Static-regridding methods are conceptually simple and robust in the sense that node movement is easy to deal with and that nodes cannot cross each other. The solution procedure can be separated into three clearly defined and independent tasks of integration, regridding and interpolation, which is advantageous in the implementation phase of the adaptive grid method. A drawback to the static-regridding method is the use of interpolation which introduces numerical errors and smooths solution profiles. Dynamic-regridding methods can compute sharp solution profiles, and may be more efficient by using fewer spatial grid points and allowing larger time steps. The coupling between the solution and the spatial variable, however, increases the size of system of equations to be solved at each time step. The moving mesh method is also prone to a phenomenon called mesh-tangling which is dealt with by the inclusion of penalty terms which govern the regularization of grid movement (Zegeling, 1991; Blom et al., 1988). Efficiency, robustness and reliability of the dynamic-regridding method depend on the successful tuning of the regularization parameters (Zegeling, 1991), and thus considerable user expertise is required. Several limitations to non-uniform grid distributions were pointed out by Thompson Chapter 3. Mesh Selection Methods 19 et al. (1985): • Nodes must concentrate, yet no region may be devoid of nodes. • Grid distribution must be sufficiently smooth and cannot be too skewed, or else the truncation error may increase. • Nodes cannot move too far or too fast, or oscillations may occur in the solution. This chapter describes the implementation of grid generating methods into the nu-merical modelling of the CBL. Two mesh selection methods are considered: a priori and static-regridding. The a priori approach is viable because the structure and development of the C B L , in terms of the mean quantities, is well understood. Static-regridding has been chosen over moving mesh methods for several reasons. Static-regridding is relatively problem independent and thus easy to implement, and computationally inexpensive. The dynamic approach, on the other hand, allots equal significance to the solution and the node positioning. Considering that second order turbulence closure models have been subject to many simplifications in order to reduce the number of prognostic equations, an important objective in this work is to avoid increasing the size of the system of equations to be solved. It should be noted that the adaptive method used in this study redistributes a fixed number of nodes in contrast to methods which add or delete nodes. Methods which redistribute nodes, but dynamically adapt the number of nodes in order to meet a user-specified tolerance (Verwer et. al., 1989) are also not considered here. The dynamic adaptation of nodes is considered to be of secondary importance, because the nodes adjust to the activity of the solution even if the number of nodes is held fixed (Verwer et al., 1989). Chapter 3. Mesh Selection Methods 20 3.1 A Priori Mesh Selection The objective is to introduce a coordinate transformation to a new independent variable ( = ((z,t) which is based on a priori knowledge of the behaviour of the solution. The transformation should be such that for the new variable, the problem can be successfully integrated on an uniform £-grid. A suitable non-uniform z-grid is then given by the inverse ( transformation z = z((,t). The coordinate transformation is based on observed characteristics of mean variables in the CBL. High spatial variability occurs in the surface and entrainment layers, while almost no variability is observed within the mixed and capping inversion layers. The desired coordinate transformation should thus introduce fine resolution in the surface and entrainment layers, and a coarser resolution in the mixed and inversion layers. In this example, the co-ordinate transformation described by Delage (1974) (log-linear transformation with a superimposed bell-shaped function which increases reso-lution around a given level) is employed. The coordinate transformation ((z) is given by CO,*) = a . , d (z — h(tY hx(z + 6) + cz + - arctan ' (3.22) e V e , where a, b, c, d, e are fixed constants. The parameter h(t) signifies the time varying loca-tion of the entrainment layer and is calculated as the height where wO is most negative. With A ( constant, the grid spacing Az is inversely proportional to d(/dz Az dz 1 d + C + (3.23) z + b (z-h(t))2 + e*\ Based on (3.23), a new grid is calculated at selected time intervals. Polynomial interpo-lation is used to transfer prognostic variables from old to new grid points. Chapter 3. Mesh Selection Methods 21 3.2 Adaptive Mesh Selection 3.2.1 Theoretical Background The concept of equidistribution implies that the error can be reduced if the grid points are distributed so that some positive weight function w{z) is equally distributed over the grid domain; ie, AziWi — constant (3.24) where Az{ = — and ze[0, H]. As in the previous section, a transformation ((z) is introduced so that the grid is uniform in (-space. By construction, let ( identify mesh points with successive integer values so that A ( = 1 and ( = 1..JV where N is the total number of nodes on the grid. Thus, Az = zrA( = zt {zc denotes dz/d() and (3.24) is now given by Zfiv = constant = C. (3.25) where the subscript i has been dropped for convenience. As discussed by Thompson et al. (1985), there are two main ways to proceed. If w is a function of (, then (3.25) is the Euler equation for the minimization of the integral which may be taken to represent the energy of a system of springs. The grid distribution resulting from equidistribution would then represent the equilibrium state of such a spring (3.26) system. The constant C is determined by integrating a variant of (3.25), = Cw~ over (e[l, N] C = H (3.27) w Chapter 3. Mesh Selection Methods 22 If the weight function, however, is a function of the grid spacing, i.e. w = w(z), then (3.25) is the Euler equation for the minimization of the integral Since ( z represents the point density, (3.28) is a measure of smoothness of the point dis-tribution. The grid distribution resulting from equidistribution represents the smoothest distribution attainable (Thompson et al., 1985). The constant C is obtained by integrat-ing a variant of (3.25), w — C(z, over ze[0,H] The above analysis shows that the constant C may be evaluated in different ways, depending on whether the weight function is associated with the mesh, i.e. £, or with the mesh point locations, i.e. z. Although both approaches are viable, the latter one, which shall be referred to as the smoothness approach, is preferable in the implementation phase of equidistribution. If the definition of C via (3.27) is used in (3.25), then the mesh adjusts to achieve a uniform value of wAz. However, this uniform value defined by (3.27) is dependent on the grid distribution and therefore an iterative approach must be taken. With the definition of (3.29), the grid adjusts to one specified uniform value since the integration occured over the physical space z. We continue then with the smoothness approach. Substituting (3.29) into (3.25) and integrating from 0 to z for z (1 to ((z) for Q, yields where z is the integral variable. Noting that C(2i) = h (3.30) may be rewritten as (3.28) (3.29) (3.30) (3.31) The values of Z{ in (3.31) represent the new grid point locations. Chapter 3. Mesh Selection Methods 23 3.2.2 B a s i c I m p l e m e n t a t i o n Equation (3.31) is solved via the variable knot spline procedure of de Boor (1974). This procedure is computationally very cheap and due to the explicit construction, the node ordering is always maintained so that the nodes do not cross each other or leave the space interval. The implementation of de Boor's technique is described in detail by Sanz-Serna and Christie (1986). The implementation is as follows (Sanz-Serna and Christie, 1986; Revilla, 1986). As-sume at time level tn a non-uniform grid zj, j — 0, . . . , J and the computed approxima-tions Uf to the true solution u(zj,tn). The solution and grid are advanced to the next time level tn+1 = tn -\- r n + i , rn+\ > 0 in two stages: 1. Use numerical scheme on grid z™ to obtain approximations U™+1 to u(z™, £ n + i ) . 2. Find a new grid z™ + 1 by equidistributing some appropriate weight function based on {J2™,I7^+1} via de Boor's technique. Compute j = 0, . . . , J , as an approx-imation to u(z" + 1, tn+i) by means of cubic polynomial interpolation of the values * j , j = o , . . . , j . It can be observed that in the implementation procedure described above, the grid lags behind the solution; that is, the solution at the (n -f l)th time level is computed on the grid at the nth time level. Verwer et al. (1989) place more emphasis on the precise grid positioning by using an intermediate step whereby the grid at the next (n + l)th time level is predicted by computing the solution at the (n + l)ih time level. The solution at the nih time level is moved to the new grid via interpolation and then the solution at (n + l)th time level is computed using the new grid. The grid is therefore ahead of the solution, but requires the system of equations to be solved twice. For our purposes, we choose the approach of Sanz-Serna and Christie (1986), since it is computationally Chapter 3. Mesh Selection Methods 24 cheaper, and it is questionable if the grid points need to be placed as accurately as Verwer et al. (1989) do. The grid distribution at t = 0 is determined by applying de Boor's algorithm to an initially uniform mesh. This process is iterated until two consecutive values of J 0 L wdz differ by less than a specified tolerance. The tolerance is set to 10 - 5 in this study. Although this iterative procedure does not guarantee convergence (Blom et al., 1988), it is assumed to provide a sufficiently accurate estimate for an initial, non-uniform mesh. 3.2.3 Weight Function From (3.24), it is apparent that the grid spacing is small where the weight function is large, and vice versa. Thus, if the weight function represents some measure of error or solution variability, then the grid spacing will be large in regions of small error or little variability, and small in regions of large error or high variability. Various forms of the weight function are discussed by Russell and Christiansen (1978) and Thompson et al. (1985). The weight function is usually taken to represent the local truncation error or some derivative of the solution. The use of formal truncation error expressions is not practical, because higher derivatives are not accurately represented through discretization and are subject to computational noise (Thompson et al., 1985). The most common approach is to equidistribute some gradients or curvature of the solution (Russell and Christiansen, 1978). In this study, the arc length AS of the norm of the mean wind velocity \\U\\ = (U2 + V2fl2 is equidistributed; that is, (AS), 2 = a{Az)j + {A\\U\\)2 = constant or ASi = yJa + (\\U\\z)?Azt '(3.32) Chapter 3. Mesh Selection Methods 25 where \\U\\Z denotes d\\U\\/dz. Comparing (3.32) to (3.24) yields the weight function The positive free parameter a can be thought of as either regulating the relative im-portance of values of z and \\U\\ (Revilla, 1986) or presenting z and \\U\\ in a suitable non-dimensional form (Sanz-Serna and Christie, 1986). \\U\\ is used as the equidistribution variable over other prognostic variables because the mean wind velocity in term of its gradient is a dominant variable in the closure parameterization and because the profiles of mean wind velocity exhibit clearly defined regions of high gradients. The mean potential temperature was excluded in the definition for the weight function, since the constant gradient of mean potential temperature in the inversion layer (lapse rate) would cause an unwanted increase in resolution in that region. (3.33) Chapter 4 Results In this chapter, the effects of various grid resolutions and distributions on model solution are evaluated by comparing the vertical profiles of mean and turbulent quantities as well as the C B L depth. Solution profiles obtained with the finest grid resolution (10m) are used as a control or reference solution. This is based on the assumption that a finer grid resolution will yield more accurate results. All profiles shown are those after five hours of simulation, unless otherwise specified. 4 . 1 Invariant Grids Before results for the mesh selection methods of Chapter 3 are presented, it is of interest to compare solution profiles for various uniform grid resolutions. It was pointed out in section 2.3, that large discrepancies in the solution profiles occur when coarse grid spacing near the surface is combined with small z\. To remove the interaction between grid spacing and z1, the only grid distributions considered here are those which are logarithmically spaced near the surface and uniformly spaced away from it. Thus, the effects of grid resolution on the entrainment layer are isolated, since the solution gradients near the surface are resolved with logarithmic spacing. Figure 4.5 shows three log-linear grids with different uniform resolution, defined as A z o o , away from the suface. The grids are computed using the co-ordinate transformation given by (3.22) with d — 0, which reduces (3.22) to a log-linear transformation. The parameters a and b are 0.5 and 0.1, respectively, and values for c are indicated in the 26 Chapter 4. Results 27 c=0.2 c=0.015 c=0.0355 J. * 0 20 40 60 80 100 120 140 AZ (m) Figure 4.5: Log-linear grid distributions for Azoo = 10, 55 and 129m corresponding to c = 0.2,0.0355 and 0.015, respectively. figure. The decrease in the last grid interval is due to a fixed upper boundary. This has negligible effects on the numerical integration since it occurs in a region of small spatial gradients. Figure 4.6 compares the vertical profiles of mean wind velocity and Reynolds stress components for Az^ = 10, 55 and 129m. The solution profiles for the two finer meshes, = 10m and Az^ = 55m, are very close. The coarser grid (Az^ = 129m) does not resolve gradients in the inversion region well and gives a lower value for the Reynolds stress component vw below the inversion, relative to the other resolutions. A slightly larger value for V in the mixed layer is also observed for the coarser mesh. The profiles for mean potential temperature and heat flux are nearly identical for all three resolutions (figure 4.7). Again, the coarser grid does not resolve the steep potential temperature or heat flux gradients in the entrainment layer. Figure 4.8 shows that the T K E profiles are similar for all three resolutions, but are overestimated by the two coarser grids in the Chapter 4. Results 28 Table 4.1: Hourly C B L height for log-linear grid distributions. Hour C B L height, Z{ (m) = 10m Azoo = 55m Azoo = 129m 1 671 688 652 2 840 851 773 3 970 960 1020 4 1089 1125 1145 5 1199 1235 1145 mixed layer as compared to the reference profile AZQQ = 10m. The entrainment heat flux w6i is approximately —0.007° K m / s for the two finer grids and —0.006° K m / s for the coarser grid. The magnitudes of these values for w$i are substantially smaller than the expected value of 0.016° K m / s , based on an entrained to surface heat flux ratio of —0.2. It is assumed that the large underestimation of wOi results from the closure parameterization, and not from the grid resolution. Table 4.1 shows the hourly C B L height for each resolution. Initially, the C B L height is 459m. The rate of change of the inversion layer, or entrainment rate, is very similiar for the finer two resolutions. The hourly averaged entrainment rate decreases from ap-proximately 210 to 110 m/hour . Discrepancies in the entrainment rate are observed for the coarser resolution wi th Az,^ = 129m, which may be due to a lack of resolution in defining Z{. For the last hour of simulation, the case Az,*, = 129m is unable to register any growth in the C B L . This result is clearly not physical, since a constant surface heat flux implies a continuously growing boundary layer. The close agreement of solution between grids Az,^ = 10m and Azoo = 55m suggests that the effects of resolution on model solution become undistinguishable beyond a cer-tain grid resolution; that is, the model solution appears to converge for ar high enough resolution. Following this argument, the solution obtained with the mesh Az,^ = 10m is Chapter 4. Results 29 Figure 4.6: Vertical profiles of mean wind velocity components (a) U and (b) V; and Reynolds stress components (c) uw and (d) vw for grid distributions shown in figure 4.5. Chapter 4. Results 30 Chapter 4. Results 31 considered to be exact for the following comparisons. 4.2 Adaptive Mesh Method For the adaptive method, the grid distribution depends on the mesh size or total number of nodes N and the parameter a, which controls the sensitivity of the grid distribution to the mean wind velocity gradient. Two grid distributions of N = 44 and 22 are considered for various a. Results obtained with these meshes are compared to the results from log-linear distributions of the same mesh size. For all graphs of vertical profiles, the results for the log-linear grid distribution with Az^ — 10m are included for comparison purposes. Figure 4.9 shows the grid distribution for N = 44, a = 0.0001 and a = 0.0005, as well as the log-linear grid of size N = 44, which corresponds to the case Az^ = 55m. Although the grid spacing for both adaptive grids is of a comparable magnitude to that of the log-linear grid of size N=44, noticable differences are found in the solution profiles. In particular, the position of the inversion base for all profiles is higher for the adaptive grids than for the log-linear grid (figures 4.10 and 4.11). This discrepancy increases with a; that is, the difference increases as the grid refinement around the entrainment layer increases. Chapter 4. Results 32 0 20 40 60 80 100 120 140 hZ (m) Figure 4.9: Grid distributions for adaptive method with N = 44, a — 0.0001 and a = 0.0005. Corresponding log-linear grid as well as reference profiles are included. - o N 8 286 287 288 289 290 291 292 -0.02 0.0 0.02 0.04 0.06 0.08 0.10 6 (K) ^0 (K m/s) Figure 4.10: Vertical profiles of (a) mean potential temperature 0; and (b) vertical heat flux w6 for grid distributions shown in figure 4.9. Chapter 4. Results 33 Figure 4.11: Vertical profiles of mean wind velocity components (a) U and (b) V; and Reynolds stress components (c) uw and (d) vw for grid distributions shown in figure 4.9. Chapter 4. Results 34 Table 4.2: Hourly C B L height for adaptive grids with N = 44. Hour C B L height, Z{ (m) a = 0.0001 a = 0.0005 1 800 730 2 965 890 3 1062 1043 4 1176 1149 5 1284 1250 Table 4.2 shows the hourly C B L height. Comparing these C B L heights to Table 4.1, it can be observed that the largest error in the entrainment rate occurs during the first hour of simulation. The grid evolution given in figure 4.12 shows that the initial non-uniform grid profile differs from the profile one hour later. The node movement necessary to adjust from the initial to the later grid profile may be a contributing factor for the accelerated entrainment rate in the first hour. a) b) — o N § o o . CM O m Z (m) 1000 / OOS OOS \ o 60 0 20 40 60 0 20 40 AZ (in) AZ (m) Figure 4.12: Grid evolution for adaptive method with N — 44, (a) a = 0.0001 and (b) a = 0.0005. Chapter 4. Results 35 c=0.2 O0.015 alpha=1 .Oe-4 alpha=5.0e-4 0 20 40 60 80 100 120 140 AZ (in) Figure 4.13: Grid distributions for adaptive method with N = 22, a = 0.0001 and a: = 0.0005. Corresponding log-linear grid as well as reference profiles are included. The grid distributions for N = 22, a = 0.0001 and a = 0.0005, as well as the log-linear grid distribution of the same size ( A ^ = 129m) are shown in figure 4.13. The height of the inversion base is overestimated for both cases of a (figure 4.14) in comparison to the reference case of AZQO = 10m. The magnitude of overestimation is comparable to the previous case of iV = 44 for the same values of a. The large discrepancies near the surface for mean wind velocity and Reynolds stress profiles (figure 4.15) are attributed to the coarse grid resolution near the surface as discussed in section 2.3. For the case a — 0.0001, oscillations develop in the mean wind velocity profiles near the top of the entrainment layer. This is not a result of the cubic interpolation scheme as the same observation is made with a quadratic scheme. The oscillation may be attributed to the large change in grid spacing for a = 0.0001. Table 4.3 shows the C B L depth at each hour of simulation for A' = 22. As for N = 44, the largest error in C B L growth occurs within the first hour of simulation. Figure 4.16 Chapter 4. Results 36 a) b) 286 287 288 289 290 291 292 -0.02 - 0.0 0.02 0.04 0.06 0.08 0.10 0 (K) ^5 (K in/s) Figure 4.14: Vertical profiles of (a) mean potential temperature 0; and (b) vertical heat flux w9 for grid distributions shown in figure 4.13. shows the grid evolution for the cases a = 0.0001 and a = 0.0005. The initial grid profiles are similar to those of the case = 44 in that very little grid refinement occurs in the entrainment region. For a = 0.0005, the grid spacing at the surface increases from 40m to 80m, thus introducing significant errors into the mean wind velocity and Reynolds stress components. To isolate the effects due to grid resolution near the entrainment layer, it would be necessary to fix the grid structure near the surface layer as in the case N = 22. The results presented above show that the solution profiles for both adaptive meshes are sensitive to the magnitude of local grid refinement around the entrainment layer, and that the discrepancies in C B L depth increase as the grid refinement increases. Sanz-Serna and Christie (1986) also point to this problem and limit the ratio of maximum to minimum grid spacing to 10. It is noted that for the adaptive meshes examined here, the ratio of maximum to minimum grid spacing is less than 3. Chapter 4. Results 37 Figure 4.15: Vertical profiles of mean wind velocity components (a) U and (b) V; and Reynolds stress components (c) uw and (d) vw for grid distributions shown in figure 4.13. Chapter 4. Results 38 Table 4.3: Hourly C B L height for adaptive grids with N — 22. Hour C B L height, Zi (m) Q = 0.0001 a = 0.0005 1 746 729 2 960 934 3 1075 1040 4 1189 1144 5 1295 1247 Figure 4.16: Grid evolution for adaptive method with N = 22, (a) a - 0.0001 and (b) a = 0.0005. Chapter 4. Results 39 4.3 A Priori Mesh Method For this method, the grid distribution is prescribed via the coordinate transformation given by (3.22). Two grids with similar ratios of maximum to minimum grid spacing but varying degree of resolution are examined: Case I: (a, 6, c, d, e) = (0.5,0.1,0.0355,300,100) Case II: (a, b, c, d, e) = (0.5,0.1,0.0150,200,100) The parameter h(t) in (3.22) defines the level around which the bell-shaped increase in resolution is centered and is calculated as the height where the turbulent heat flux is most negative. The results are presented for various regridding frequencies and interpolation schemes. The regridding frequencies are frg = 1 and 120 which correspond to regridding at every time step and every 120ih time step, respectively. Polynomial interpolation schemes of second and third order are tested. Results are compared to profiles obtained with corresponding log-linear grids of identical Az^, rather than grids of the same mesh size. In addition, results for the log-linear grid distribution with AZoo — 10m are included for reference. The grid distributions after 5 hours are shown in figure 4.17 for Case I. It is noted that the grid profiles for frg = 120 overlap. Figure 4.18 shows the corresponding profiles for mean wind velocity and Reynold stress components. Oscillations in the mean wind velocity components are observed near the inversion base for both simulations using cu-bic polynomial interpolation. The oscillations are much more pronounced for the higher regridding frequency (fTg — 1). Quadratic polynomial interpolation removes the oscilla-tions, but smooths out the profiles of the mean wind velocity components at the bottom of the inversion layer for the case fTg — 1. The large deviations from the reference profiles observed for the Reynolds stress components with quadratic interpolation and fTg = 1 as Chapter 4. Results 40 control log-linear ignd=1 ,cu igrid=120,cu igrid=1 ,qu igrid=120,qu 0 20 40 60 80 100 120 140 AZ (m) Figure 4.17: Grid distributions for a priori method, Case I. fTg is the regridding frequency; cubic interpolation is indicated by cu and quadratic interpolation by qu. Corresponding log-linear grid as well as reference profiles are included. well as the smaller irregularities in the Reynolds stress and heat flux profiles (figure 4.19b) for all other cases are believed to be artifacts of the interpolation. Figure 4.19a shows that the mean potential temperature profiles are almost identical for all cases, but that 0 is underestimated in the mixed layer and the inversion layer is substantially higher as compared to the reference case. Table 4.4 shows that the entrainment rate is affected by the regridding frequency, but not by the interpolation scheme. Again, the largest difference in the entrainment rate as compared to the reference solution occurs in the first hour. In contrast to the adaptive method, however, the magnitude of grid refinement around the entrainment layer is constant during the entire simulation and does therefore not explain the accelerated entrainment rate. The grid distributions for Case II after 5 hours and the corresponding log-linear Chapter 4. Results •41 Figure 4.18: Vertical profiles of mean wind velocity components (a) U and (b) V; and Reynolds stress components (c) uw and (d) vw for grid distributions shown in figure 4.17. Chapter 4. Results 42 ^ ) b ) 286 287 288 289 290 291 292 -0.02 0.0 0.02 0.04 0.06 0.08 0.10 0 ( K ) ( K m / s ) Figure 4.19: Vertical profiles of (a) mean potential temperature 0; and (b) vertical heat flux w6 for grid distributions shown in figure 4.17. Table 4.4: Hourly C B L height for a priori grids, Case 1. Hour C B L height, Zi (m) frg — 1 frg = 120 cubic quadratic cubic quadratic 1 750 750 741 741 2 919 919 958 958 3 1074 1074 1116 1116 4 1218 1183 ' 1233 1233 5 1327 1317 1343 1343 Chapter 4. Results 43 grid with Az^ — 129m are shown in figure 4.20. Similar to Case I, the grid profiles for frg = 120 are identical, but a much larger spread in the C B L height is observed for frg = 1. Figure 4.21 further shows the large effects of regridding frequency and interpolation scheme on the solution. The mean wind velocity profiles obtained with the higher regridding frequencey (frg — 1) are nonphysical in that the inversion base is not defined for the quadratic case. For the cubic interpolation, no entrainment layer is visible for the mean wind velocity component U while V exhibits an undulation which extends to the upper boundary. The oscillations in the turbulent fluxes for the case of cubic interpolation and frg = 1 are likely caused by the effect of the upper boundary condition on V propagating down into the boundary layer. For the lower regridding frequency (frg — 120), virtually no differences are found in the profiles of both mean wind velocity and Reynolds stress components for the various interpolation schemes. For all profiles of mean wind velocity as well as mean potential temperature (figure 4.22), the inversion height is overestimated by approximately 200m as compared to the reference solution. The results presented above indicate that for the a priori method the node placement affects the accuracy of the interpolation scheme and that the error introduced into the solution via the interpolation is cumulative. Chapter 4. Results 44 0 20 40 60 80 100 120 140 AZ (m) Figure 4.20: Grid distributions for a -priori method, Case II. frg is the regridding fre-quency; cubic interpolation is indicated by cu and quadratic interpolation by qu. Corre-sponding log-linear grid as well as reference profiles are included. Chapter 4. Results 45 a) b) -0.4 -0.3 -0.2 -0.1 0.0 0.1 -0.15 -0.10 -0.05 0.0 0.05 0.10 uw (m2/s2) vw (m'2/s2) Figure 4.21: Vertical profiles of mean wind velocity components (a) U and (b) V; and Reynolds stress components (c) uw and (d) vw for grid distributions shown in figure 4.20. Chapter 4. Results 46 Figure 4.22: Vertical profiles of (a) mean potential temperature 0; and (b) vertical heat flux w6 for grid distributions shown in figure 4.20. Chapter 5 Summary of Results and Conclusions This thesis presents the first implementation of time-adaptive grids in higher order turbu-lence closure modelling of the CBL. Two numerically inexpensive grid generating meth-ods, a priori and adaptive, are tested and compared against invariant log-linear grids of the same mesh size as well the reference mesh, given by Az^ = 10m. The results are summarized below: • A comparison between log-linear grids of various mesh sizes shows almost identical results for the grid Azoc — 55m and the reference mesh. Larger discrepancies, especially in the entrainment layer, are found for the case Az^ = 129m. • All grid generating methods tested in this study overestimate the height of the inversion layer. • The error in the hourly entrainment rate is largest for the first hour of simulation for all non-uniform grids. • For both the a priori and the adaptive mesh methods, the discrepancies in model solution as compared to the reference solution decrease as the mesh size increases. • For both methods, the solution error increases as the grid refinement around the entrainment layer increases. • Solutions obtained with the adaptive mesh method do not depend on the interpo-lation scheme. 47 Chapter 5. Summary of Results and Conclusions 48 • The a priori mesh method is sensitive to the regridding frequency as well as the interpolation scheme. The sensitivity to the interpolation scheme increases as the regridding frequency is increased. No differences in the solutions are observed for the lower regridding frequency. The discrepancies caused by frequent regridding are larger for smaller mesh sizes. It is concluded that the non-uniform grid distributions investigated in this thesis presently do not provide any advantages over uniform meshes of the same size due to the inability to correctly predict the height of the inversion layer. The applicability of these methods in C B L modelling is further disputed by the sensitivity of the solutions to mesh size and local grid refinement. The close agreement in results between the log-linear grids Az^ = 10m and Az^ = 55m, and to some extend Azoo = 129m, suggests that grid resolution does not greatly influence the solutions obtained with the simple Level 2.5 turbulence closure model. The importance of grid resolution on solution accuracy may depend on the complexity of the physical model; that is, results obtained with higher level turbulence closure models may be more sensitive to grid resolution. Possible extensions to this work include further examinations of the effects of inter-polation on the a priori mesh method as well as the effects of the initialization procedure of prognostic variables on both grid generating methods. Bibliography Andre, J . , G. de Moor, P. Lacarrere, G. Therry and R. du Vachat, 1978: Modeling the 24-hour evolution of the mean and turbulent structures of the planetary boundary layer. J. Atmos. Sci., 35, 1861-1883. Andren, A. , 1990: Evaluation of a turbulence closure scheme suitable for air-pollution applications. Journal of Applied Meteorology, 29, 224-239. Ascher, U. M . , R. M . Mattheij and R. D. Russell, 1988: Numerical Solution of Boundary Value Problems for Ordinary Differential Equations. Prentice-Hall series in computa-tional mathematics. Prentice Hall. Blackadar, A. , 1962: The vertical distribution of wind and turbulent exchange in neutral atmosphere. J. Geophy. Res., 67, 3095-3102. Blom, J . , J . Sanz-Serna and J. Verwer, 1988: On simple moving grid methods for one-dimensional evolutionary partial differential equations. Journal of Computational Physics, 74, 191-213. Cai, X. , 1993: Large Eddy Simulation of the Upper Atmospheric Surface Layer. Ph.D. thesis, University of British Columbia. Caughey, S., 1982: Observed characteristics of the atmospheric boundary layer, in F. Nieuwstadt and H. van Dop, editors, Atmospheric Turbulence and Air Pollution Modelling, Atmospheric Science Library, pp. 107-158. D. Reidel Publishing Company. 49 Bibliography 50 de Boor, C , 1974: Good approximation by splines with variable knots. II. in G. Watson, editor, Conference on the Numerical Solution of Differential Equations. Dundee 1973, Vol. 363 of Lecture Notes in Mathematics, pp. 12-20. Springer Verlag. Deardorff, J. W., 1972: Numerical investigation of neutral and unstable planetary bound-ary layers. J. Atmos. Sci., 29, 91-115. Delage, Y. , 1974: A numerical study of the nocturnal atmospheric boundary layer. Quart. J. R. Met. Soc, 100, 351-364. Driedonks, A. and H . Tennekes, 1984: Entrainment effects in the well-mixed atmospheric boundary layer. Boundary-Layer Meteorology, 30, 75-105. Garratt, J . , 1992: The atmospheric boundary layer. Cambridge Atmospheric and Space Science Series. Cambridge University Press. Helfand, H . and J. Labraga, 1988: Design of a nonsingular level 2.5 second-order closure model for the prediction of atmospheric turbulence. J. Atmos. Sci., 45(2), 113-132. Kolmogoroff, A. , 1941: The local structure of turbulence in incompressible viscous fluid for very large reynolds number. Dokl. Akad. Nauk. SSSR., 30, 301. Lumley, J . and B. Khajeh-Nouri, 1974: Computational modeling of turbulent transport. Advances in Geophysics, 1 8 A , 169-192. Mellor, G. L. and T. Yamada, 1974: A hierarchy of turbulence closure models for plane-tary boundary layers. J. Atmos. Sci., 31, 1791-1806. Mellor, G. L. and T. Yamada, 1982: Development of a turbulence closure model for geophysical fluid problems. Reviews of Geophysics and Space Physics, 20(4), 851-875. Bibliography 51 Miyakoda, K. and J. Sirutis, 1977: Comparative integrations of global models with var-ious parameterized processes of subgrid scale vertical transports. Beitr. Phys.Atmos., 50, 445-485. Paulson, C , 1970: The mathematical representation of wind speed and temperature in the unstable atmospheric surface layer. J. Appl. Meteor., 9, 857-861. Revilla, M . , 1986: Simple time and space adaptation in one-dimensional evolutionary partial differential equations. International Journal for Numerical Methods in Engi-neering, 23, 2263-2275. Rotta, J . , 1951: Statistische Theorie Nichthomogener Turbulenz. Z. Phys., 129, 547-572. Russell, R. and J . Christiansen, 1978: Adaptive mesh selection strategies for solving boundary value problems. SIAM J. Numer. Anal, 15(1), 59-80. Sanz-Serna, J. and I. Christie, 1986: A simple adaptive technique for nonlinear wave problems. Journal of Computational Physics, 67, 348-360. Stull, R. B., 1988: An Introduction to Boundary Layer Meteorology. Atmospheric Science Library. Kluwer Academic Publishers. Taylor, P. and Y . Delage, 1971: A note on finite-difference schemes for the surface and planetary boundary layers. Boundary-Layer Meteorology, 2, 108-121. Tennekes, H. and A. Driedonks, 1981: Basic entrainment equations for the atmospheric boundary layer. Boundary-Layer Meteorology, 20, 515-531. Thompson, J . F. , 1985: A survey of dynamically-adaptive grids in the numerical solution of partial differential equations. Applied Numerical Mathematics, 1, 3-27. Bibliography 52 Thompson, J. F. , Z. Warsi and C. W. Mastin, 1985: Numerical Grid Generation: Foun-dations and Applications. North-Holland. Verwer, J . , J. Blom and J. Sanz-Serna, 1989: An adaptive moving grid method for one-dimensional systems of partial differential equations. Journal of Computational Physics, 82, 454-486. Yamada, T., 1977: A numerical experiment on pollutant dispersion in a horizontally-homogeneous atmospheric boundary layer. Atmos. Environ., 11, 1015-1024. Yamada, T. , 1983: Simulations of nocturnal drainage flows by a q2l turbulence closure model. J. Atmos. Sci., 40, 91-106. Yamada, T. and G. Mellor, 1975: A simulation of the wangara atmospheric boundary layer data. J. Atmos. Sci., 32, 2309-2329. Zegeling, P., 1991: Moving-Grid Methods for Time-Dependent Partial Differential Equa-tions. CWI Tract 94. Centrum voor Wiskunde en Informatica. Zeman, 0., 1981: Progress in the modeling of planetary boudary layers. Ann. Rev. Fluid Mech., 13, 253-272. Zeman, 0. and J. L. Lumley, 1976: Modeling buoyancy driven mixed layer. J. Atmos. Sci., 33, 1974-1988.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A convective boundary layer model with time-adaptive...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A convective boundary layer model with time-adaptive grids Rucker, Magdalena 1995
pdf
Page Metadata
Item Metadata
Title | A convective boundary layer model with time-adaptive grids |
Creator |
Rucker, Magdalena |
Date Issued | 1995 |
Description | This thesis presents the first implementation of time-adaptive non-uniform grid distributions in higher order turbulence closure modelling of the convective boundary layer. The aim was to introduce a vertical grid that resolves the surface and entrainment layers, but uses coarser resolution elsewhere. Two types of grid generating methods were examined with the Level 2.5 turbulence closure model (Mellor and Yamada, 1982; Yamada, 1983). The solution profiles of mean and turbulent quantities were compared to model results obtained with stationary grid distributions of the same mesh size as well as a reference solution obtained with a grid resolution of 10m. |
Extent | 2353931 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-01-10 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0079747 |
URI | http://hdl.handle.net/2429/3526 |
Degree |
Master of Science - MSc |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1995-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1995-0074.pdf [ 2.24MB ]
- Metadata
- JSON: 831-1.0079747.json
- JSON-LD: 831-1.0079747-ld.json
- RDF/XML (Pretty): 831-1.0079747-rdf.xml
- RDF/JSON: 831-1.0079747-rdf.json
- Turtle: 831-1.0079747-turtle.txt
- N-Triples: 831-1.0079747-rdf-ntriples.txt
- Original Record: 831-1.0079747-source.json
- Full Text
- 831-1.0079747-fulltext.txt
- Citation
- 831-1.0079747.ris
Full Text
Cite
Citation Scheme:
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>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0079747/manifest