THE DEVELOPMENT AND APPLICATION OF A ONE-LEVEL MESOSCALE ATMOSPHERIC MODEL OVER THE LOWER FRASER VALLEY OF BRITISH COLUMBIA by KEITH WILLARD AYOTTE B.Sc. (Hons.) The University of British Columbia, 1984 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE STUDIES (Department of Geography) We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA August 6,1986 © Keith Willard Ayotte In presenting this thesis in p a r t i a l f u l f i l m e n t of the requirements for an advanced degree at the University of B r i t i s h Columbia, I agree that the Library s h a l l make i t 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 i s understood that copying or publication of t h i s thesis for f i n a n c i a l gain s h a l l not be allowed without my written permission. Department of Geogr^ ph-* The- University of B r i t i s h Columbia 1956 Main Mall Vancouver, Canada V6T 1Y3 ^ ^ 26/August/1986 Date -6 (3/81) - i i -ABSTRACT This study involves the use of a one-level mesoscale numerical model to simulate thermally driven, topographically chanelled atmospheric flow over the Lower Fraser Valley of British Columbia. The model used is a modified version of a model developed by Mass Dempsey at the University of Washington. The model was modified to include friction parameterization that was spatially variable and stability dependent. The diabatic surface heating parameterization was also modified to account for spatial variability of surface thermal characteristics. Model runs of July 20,1985 and August 23,1985 are presented in this study. The model is validated using data from a measurement network over the model domain. Statistics of the validation on the two days modeled are presented. It was found that the modified model is quite sensitive to the surface thermal characteristics specified but relatively insensitive to spatial variability of surface roughness. The model is limited in its ability to resolve effects of small topographic features in high Froude number situations but the large scale flow around topographic features is well modeled. Synoptic scale thermally induced pressure gradients are found to affect ability of the model to simulate near surface flow within the domain . It is suggested that the model be nested in a lower resolution grid in order to account for this forcing. It was concluded that this modification, coupled with more accurate specification of the surface characteristics over the model domain would enhance the agreement between modeled and observed fields. T A B L E OF CONTENTS Abstract Table of Contents L i s t of Tables L i s t of Figures Acknowledgements 1. Introduction 2. The Model 2.1 Model Formulat ion 2.2 Model Equation Solutions 2.2.1 Init ialization 2.2.2 Evolut ion Through Time of Modeled Fields 2.3 Model Outputs 3. Surface Parameterization 3.1 Diabatic Heating/Cooling Parameter izat ion 3.2 Surface Frict ion Parameterizat ion 4. Exploratory Modeling 4.1 Testing of Surface Frict ion Parameter izat ion 4.2 Test R u n Over a Simple Topographic Feature 4.3 Model ing Over Idealized Topography 4.3.1 Edge Fi l ter Effects 5. Modeling Over The Real Domain Surface 5.1 U p p e r Boundary Conditions 5.1.1 J u l y 20,1985 5.1.2 A u g u s t 23,1985 -iv-5.2 Statistics used in Model Validation 68 5.3 Modeled and Observed Fields and Validation Data 69 5.4 Discussion of Modeling of July 20,1985 70 5.5 Discussion of Modeling of August 23,1985 115 5.6 General Observations 156 6. Conclusions and Recommendations for Further Work 161 Bibliography 165 Appendix I: Sigma-Coordinate Transformation 173 Appendix II: Derivation of the Tendency Equations 177 Appendix III: Collection of Validation Data 182 Appendix IV: Radiation Input and Slope Correction 188 Appendix V: Calculation of Temperature Tendency from Surface Sensible Heat Flux 204 Appendix VI: Statistical Indices 209 Appendix VII: Model Run Parameters 210 Appendix VIII List of Symbols 211 - V-LIST OF TABLES page Table 4.1.1 Friction Parameterization Comparison 48 Table a3.1 List of Anemometer Locations 180 - v i -LIST OF FIGURES page 1.1 The Lower Fraser Valley Region 2 2.1.1 Vertical Temperature Structure Parameterization 14 3.1 Scales of Atmospheric Motion 24 3.1.1 Bowen Ratio Over the Modeled Domain 28 3.1.2 Modeled Energy Budget Terms, July 20,1985 29 3.1.3 Modeled Energy Budget Terms, August 23,1985 30 3.1.4 Measured Energy Budget terms 31 3.1.5 Sensitivity of Surface Pressure Change to Height of Topographic Influence 32 3.2.1 Roughness Length Over the Modeled Domain 38 3.2.2 U * Stability Correction - Exact vs Approximated 42 3.2.3 U * Stability Correction vs Time -July 20,1985 43 3.2.4 U * Stability Correction vs Time - August 23,1985 44 4.1.1 Frictional Turning - Modeled vs Calculated 49 4.2.1 Sinusoidal Hill Topography 51 4.2.2 Near Surface Flow Around a Sinusoidal Hill 52 4.2.3 Near Surface Temperature Pattern Over a Sinusoidal Hill 53 4.3.1 Real and Idealized Topography 57 4.3.2a Flow Over Idealized Topography - 06:00 58 4.3.2 Flow Over Idealized Topography - 12:00 59 5.1.1 850mb Chart- 1200Z July 20,1985 66 5.1.2 850mb Chart - 1200Z August 23,1985 67 5.4.1a-24a Modeled Near Surface Winds - July 20,1985 78-113 5.4.1b-24b Observed Near Surface Winds - July 20,1985 78-113 - v i i -5.4.1c-24c Wind Velocity Residual Plots - July 20,1985 79-112 5.4.25 Indices of Agreement vs Time - July 20,1985 114 5.5.1a-24a Modeled Near Surface Winds - August 23,1985 119-154 5.5.1b-24b Observed Near Surface Winds - August 23,1985 119-154 5.5.1c-24c Wind Velocity Residual Plots - August 23,1985 120-153 5.5.25 Indices of Agreement vs Time - August 23,1985 155 5.6.1 Coefficient of Variation vs Wind Speed July 20,1985 159 5.6.2 Coefficient of Variation vs Wind Speed August 23,1985 160 a l . l Constant Pressure Surfaces Over Topography 175 a 1.2 Constant Sigma Surfaces Over Topography 176 a3.1 Validation Data Collection Sites 181 a4.1 Idso and Jackson (1969) Net Longwave Radiation Parameterization 189 4.2 Direct Normal Solar Radiation - July 20,1985 193 a4.3 Direct Solar Radiation - July 20,1985 194 a4.4 Global Solar Radiation - July 20,1985 195 a4.5 Diffuse Solar Radiation - July 20,1985 196 a4.6 Anisotropy Index - July 20,1985 197 a4.7 Direct Normal Solar Radiation August 23,1985 198 a4.8 Direct Solar Radiation - August 23,1985 199 a4.9 Global Solar Radiation - August 23,1985 200 a4.10 Diffuse Solar Radiation - August 23,1985 201 a4.11 Anisotropy Index - August 23,1985 202 a4.12 Topography Defined by Grid Node Elevation 203 a5.1 Elevation vs Modeled Temperature Tendency 208 - v i i i -ACKNOWLEDGEMENTS The debt owed those who supported this research goes far beyond what can be expressed in a few short pages. I would like thank m y supervisor D r . Douw G . Steyn for his technical and organizational support throughout this research. In addition, D r . Stejm has provided much advice and counselling on a number of academic points s tart ing a number of years prior to m y entr j ' into graduate school. F o r this I a m deeply grateful . I would also like to thank D r . T imothy R. Oke and D r . John E . H a y for their patient and insightful guidance on the research carried out for this thesis. M a n y thanks are also due D r . D a v i d D e m p s y and D r . Cli f ford M a s s for the use of their model in this s tudy. Without their contribution the project could not have gone ahead. M y thanks go to D r . Ian M c K e n d r y for his help in the preparation of model validation data, Scott Robeson for his help with the statistical aspects of this stud}' and Stephen R a m s a y for laborator j ' assistance in the earliest stages of the Tethersonde program. Assistance in the collection of data was provided by M a x B r o w n . The research was supported by grants, subventions and contracts to D r . D . G . Steyn by the N a t u r a l Science and Engineering Research Council of Canada, the Atmospheric Envi ronment Service of C a n a d a and the B r i t i s h Columbia M i n i s t r y of Environment . M y personal funding was provided as Research Assistantships from the above sources, and Teaching Assis tantships f rom the Department of Geography at the U n i v e r s i t y of Br i t i sh Columbia . The operation of the Tethersonde was made possible by the courteous and efficient efforts of M r . Gale Winterburn of Transport Canada. Permission to locate the anemometry network instruments on commercial , industr ia l and institutional locations was graciously given by: M r . R. M a y of M a y Brothers F a r m s , M r . J . Lucket t of the Atmospheric Envi ronment Service, L t - C d r . M . E . Stanley of C F S Aldergrove, M r . R. Ruock of Ex t raord ina i r , M r . H . Lee of the M u n i c i p a l i t y of Richmond, M r . I. Vinegar of M a s t e r w a s h - i x -Products and M r . R . G . Gascoyne of the C i t y of White Rock. Addit ional data were provided by M r . J . Lucket t of the Atmospheric E n v i r o n m e n t Service, D r . R . C . Bennett of the B r i t i s h Columbia M i n i s t r y of Environment and M r . M . Mennel l of the Greater Vancouver Regional Distr ic t . Assistance f rom Tom Nicol and J i m Glosl i of the Computing Centre at the U n i v e r s i t y of B r i t i s h Columbia in us ing the array processor is greatly appreciated. L a s t but far f rom least, I owe a great deal to Jane Roots for her m a n y skil ls used in the preparation of this document and letting sanity prevai l around the home. -1-1. I N T R O D U C T I O N The main objective of this study is the investigation of the use of a one-level mesoscale numerical model to simulate windfields over topography having complex relief and surface characteristics. Within this framework there are a number of aspects of the modeling addressed that concern the model itself and the assumptions under which it operates. Among these aspects are the model formulation and the representation of the real atmosphere provided by that formulation as well as the manner in which the model handles spatially varied surface characteristics. It is also the intent of the author to establish a means by which the model can be objectively validated. This is done for two reasons. The model results must be validated in order to assess the performance of the model and to aid in the development of the model parameterizations. The other equally important reason is to establish a basis for comparison of the work presented here to any further work done with this model or with other attempts to model the same fields using different methods. The work undertaken here was carried out at and around the University of British Columbia in Vancouver. This study is part of a long term project undertaken by a number of individuals under the supervision of Dr. D.G. Steyn in the Department of Geography at the University of British Columbia. The objective of this encompassing project is the investigation of the mesoscale flow patterns over the lower Fraser Valley of British Columbia as well as the effects of these flows and circulations on the air quality in the area. Vancouver, British Columbia is a city of approximately 1.5 million people including the surrounding municipalities. The city is built on the large deltaic deposit which has the Strait of Georgia as the western boundary and extends some 120 kilometres to the east-north-east where the valley narrows to enter the Fraser canyon, (see Figure 1.1) To the north the flat valley bottom is bordered by mountains rising to just over 1000 metres. In plan, the valley has a general triangular shape with mountains on the south-eastern edge - 3 -similar in size and shape to those of the northern boundary of the valley. This rather picturesque setting is the host to many interesting atmospheric phenomena. The Fraser Valley lies at approximately 50 degrees North latitude and, in the winter months, is affected by an almost continual sequence of mid-latitude disturbances. These disturbances from the north-west and south-west bring with them a large amount of precipitation usually in the form of rain. The mean annual rainfall varies over the basin from a low in the south-western corner of 1000 mm, increasing to the north to over three times this amount as the topography rises. The precipitation at these higher elevations is largely in the form of snow. (Hay and Oke,1976) This orographic enhancement by the mountains on the northern boundary of the valley is only one of the many phenomena induced by the unique topographic setting of the lower Fraser Valley. There exists a very sharp contrast between the ocean surface to the west and the flat bottom of the valley in terms of thermal and roughness characteristics. Though not as evident as the sea-surface/valley bottom contrast, similar variations also exist within the valley itself. These variations and sharp contrasts, acting with the local and regional topography, can produce some very interesting atmospheric phenomena. For example, in the summer months when the synoptic scale disturbances are steered to the north by a blocking high, the synoptic pressure gradient over the lower Fraser Valley is very slack. The high solar loading along with the nearly non-existent synoptic pressure gradient provides the setting for the formation of thermally forced topographically channeled mesoscale flows and circulations. One of the most predominant of the thermally driven atmospheric features is the land/sea breeze circulation. Given the high mountain boundaries to the north and south-east and the circulatory nature of the land/sea breeze flow, the ventilation of atmospheric pollutants in the lower Fraser Valley can be very poor during a prolonged period of high surface pressure and the associated slack synoptic pressure gradient. This problem is exacerbated by the limited vertical mixing accompanying these circulations. This has some -4-rather startling implications in terms of air quality within the basin. The build-up of oxidants due to the reaction of precursor hydrocarbons and oxides of nitrogen in the presence of high levels of ultra-violet light can present health hazards as well as detract from the natural beauty of the setting. (Concord Scientific, 1982,1985; Wilson, Mills and Wituschek,1984) From a more broad, somewhat more academic viewpoint, flows and circulations of this nature are of great interest. These phenomena are, in part, the response of the atmosphere to spatially varying surface characteristics and complex topography. For example Bornstein and Johnson (1977), in an examination of urban-rural wind velocity differences, present data which clearly shows the effects of varying surface roughness and thermal characteristics. Hjelmfelt (1982) examines urban versus non-urban effects to mesoscale boundary layer airflow and vertical air motion. Hjelmfelt suggests that in certain cases non-urban effects interact with urban-induced vertical air motion and that topograph}' is the dominant factor in producing the major non-urban effects. Hjelmfelt adds to this: "For different cities with differing urban characteristics and differing geographical settings, their net effects in combination may be different". Channeling of flow around complex topography has been examined in studies carried out within domains that are near or include the lower Fraser Valley (Overland and Walter, 1981,1983; Schoenberg,1983; Mass, 1981,1982; Danard 1977). All of these studies are on a much larger scale than that of the lower Fraser Valley and its immediate surroundings and none address the effects of spatially varying surface characteristics. Vasanji and Gartshore (1978) acknowledge the existence of "local terrain effects" in a study carried out in the western-most section of the lower Fraser Valley. In this study no attempt was made to establish a causal relationship between these "local terrain effects" and observed spatial variability in wind data. If Hjelmfelt's (1982) suggestion is correct and the interactions between varying surface characteristics (urban/non-urban) and complex topography may -5-be different in different settings, there is yet a great deal to be learned about these effects in the setting chosen for this study. To the extent that a model of the atmosphere can represent the atmosphere's response to variations in surface characteristics and complex topography, a model can be used to help understand these effects. If these features can be modeled and validated in an objective manner, then the model may serve as a tool to test the response of the atmosphere to inhomogeneities in roughness and thermal characteristics. The first step in the modeling process, and the one with which the major part of this thesis is concerned, is the development of a model suitable for the task of studying these flows. There are a number of constraints on the model invoked by physical limitations and, to a lesser extent, computing resource limitations. At the same time the model must produce reasonably accurate time dependent near-surface windfields and it must do so at a low cost in terms of computer resources and input data required. Two types of models present themselves as possible tools for diagnosing time dependent near-surface windfields; the fully three-dimensional model and the one-level model that uses a parameterization of the vertical structure of the atmosphere. The fully three-dimensional model, because it explicitly calculates upper level fields, may yield a more accurate near-surface windfield. Though the upper level fields are of interest in an indirect way, their calculation is not one of the main objectives of this investigation. If only near-surface windfields are required, the use of a fully three-dimensional model seems to be less than efficient in achieving this goal. A one-level model has the advantage of not having to explicitly compute upper level fields. This in turn means that the model uses much less computer time to produce the near-surface windfields desired. In the choice of a one-level model, a compromise is accepted. The one-level model uses less computer time to run but may be less accurate than its three-dimensional counterpart. - 6 -One-level models have been used with success, to examine near-surface flow over complex terrain. (Danard, 1971,1977; Danard and Thompson, 1983; Lavoie, 1972,1974; Mass, 1981,1982; Mass and Dempsey, 1985a, 1985b) These models are quite economical to operate, and within certain limitations, do a good job of modeling near-surface windfields. Some of the limitations to these models are found in the parameterization of the vertical structure and as such will be discussed in the section that deals with this parameterization. Other limitations are of a more general nature and involve the extent and the manner in which the model represents the physics of the atmosphere. In order to better understand these limitations and put them into context, the development of schemes for diagnosing surface windfields can be examined. One of the main points of this examination is to see how the more complex physically based models have grown from the basic interpolation schemes used in the first attempts to diagnose surface flow. 1.1 The Development of Objective Analyses Since as early as 1911, the method to account for spatial variability in scalar and vector fields by the use of a simple mathematical model has been known (Theissen,1911). It was realized that the interpolated values-of precipitation needed to be weighted before they were averaged to avoid biasing due to irregular spacing of observation points. Early efforts at objective analysis involve many, often complicated weighting schemes. Goodin, McRae and Seinfeld (1979) present a survey of the early methods of interpolation of sparse data, in which the evolution of weighted, least squares polynomial and optimum interpolation are presented. With the increase in the speed and availability of computer resources, the ability of the analyst to find solutions to problems of objective analysis has grown and taken a different tangent. Where in the past the solutions were arrived at via simple weighting schemes, they have more recently been found by more and more complex schemes. The evolution - 7 -has turned in the direction of physically based analyses which make use of the actual modeling of fields through the numerical simulation of their behavior. Implicit in these simulations is the deterministic nature of the phenomena of interest and the adequate understanding of the immanent processes involved in their occurrence. 1.2 Models Used in Physically Based Analyses Endlich (1967) developed a method for objective analysis that involves an iterative method to alter the kinematic properties of the windfield. In this method, iterative adjustment of the windfield is done subject to constraints on vorticity and divergence. A finite-difference scheme is used to calculate the divergence and vorticity. Endlich's method ultimately uses "meteorological knowledge" to guide the final results, making the model somewhat subjective but still physically based. Anderson (1971), Dickerson (1978) and Sherman (1978) make use of models based on conservation of mass in two or three dimensions to calculate mass-consistent windfields. Anderson (1971) uses least-squares techniques over an interpolated windfield to yield a mass-consistent windfield whereas Dickerson (1978) and Sherman (1978) apply variational methods to create a mass-consistent three-dimensional windfield model. This variational technique was developed by Sasaki (1970) and defines a functional whose external solution minimizes the difference between observed and calculated values subject to constraints. The derived Euler-Lagrange equations are then solved by finite-difference methods to minimize the functional. This method uses strong and weak constraints; weak constraints being applied approximately and strong constraints being applied exactly.These models have been used in a number of locations with a high degree of success and appear to be practical tools within an air quality framework. Goodin, McRae and Seinfeld (1979) employ a method of analysis using a three-dimensional windfield model. This model uses terrain-following coordinates and variable vertical grid spacing. It is initialized by interpolating the surface and upper wind -8-measurements. A local terrain adjustment technique involving the solution of a Poisson equation is used to establish the horizontal component of the near-surface windfield. Vertical velocities are derived from successive solutions of the continuity equation, followed by an iterative procedure which reduces anomalous divergence in the complete windfield. To this point in the development presented here, the analysis schemes have been interpolative in the strictest sense. This means that to compute a windfield over a domain consisting of a regularly spaced grid, there must first exist scattered windfield data over the domain. Another method for the interpolation of meteorological parameters at the surface is the integration of the equations that govern the state of the atmosphere. Specifically, these are the solutions to initial or boundary value problems that are posed in such a way as to include the relevant physics of the phenomena of interest. Given the complexity and the inherent non-linearity of the equations involved, use of numerical methods to integrate these equations is necessary. The exceptions to this are linearized, very much simplified sets of equations which can be solved analytically. These are seldom of use in the modeling of real windfields. (Pielke,1984) The form and complexity of the equations as well as the resolution at which the equations are applied is dependent upon the phenomena being examined and the computing resources available. There are in existence a number of three-dimensional mesoscale models that can simulate such phenomena as nocturnal low level jets, squall lines, mountain valley circulations, land/sea breezes and urban heat island flows (Keyser and Anthes,1977; Pielke,1984). The main drawback to the implementation of these models, at least in an operational mode, is the prohibitive amount of computing power needed. The high cost of running these models necessitates the striking of a balance between simplicity and the extent and manner in which the model represents the physics of the atmosphere. Through scale analysis, a number of simplifying assumptions can be made to permit the solutions of the governing equations to be found in a much less difficult and more economical fashion. For example the "shallow water" equations are arrived at by -9-assuming that the atmosphere is incompressible. This assumption is valid if the depth of the circulation is much less than the scale of the deep atmosphere (Pielke, 1984; Haltiner and Williams, 1980). This modified version of the equations is much less complex than the more general equation in that a number of terms are neglected, but also has a more restricted application. The above is only one of a number of assumptions that can be made to simplify a mesoscale model (Anderson, Tannehill and Pletcher,1984). Each one of these assumptions will in turn restrict the use of the model by ignoring aspects of the physics of the atmosphere (Wippermann,1981). It is here that part of the balance between simplicity and physics must be struck. There are also a number of constraints which guide the choice of model formulation that are for the most part numerical in origin. As an example, the transfer of energy from low wavenumbers to high wavenumbers, when modeled on a grid, can lead to non-linear instability (Pielke, 1984). This instability can be a result of a mismatch of scales between phenomena in the real atmosphere and the resolution and numerical formulation of the model. This mismatch of scales can most easily be seen in the the cascade of energy from low to high wavenumbers. Energy entering the atmosphere, does so at a wide spectrum of length and time scales. Kinetic energy can exist in the form of perturbations in the atmosphere ranging from the very low wavenumbers (planetary waves) to the very high wavenumbers (small scale turbulence). The energy at low wavenumbers tends to cascade down to higher wavenumbers. This cascade of energy continues down through motions of the length and time scales of the inertial sub-range of turbulence, to finally be removed from the atmosphere by viscous dissipation (Beer, 1974). One of the mechanisms that facilitates this cascade of energy is non-linear wave-wave interaction (LeBlond and Mysak,1978). The result of two waves interacting can be waves, of significantly higher and lower wavenumbers. The limitations of a numerical model are such that it cannot resolve - 1 0 -wavelengths smaller than the length of two model grid spaces. If the high wavenumber waves created by the wave-wave interaction in the real atmosphere are beyond this limit then thej' cannot be represented by the model. The model's inability to represent this part of the cascade may result in aliasing, where the energy from the high wavenumbers folds back into the lower wavenumbers. This can lead to inaccuracy in the modeled results or even the exponential growth of erroneous modeled fields. This is clearly undesirable. To prevent this aliasing, smoothing filters may be applied to the modeled fields or the model inputs in order to prevent this numerical instability (Pielke,1984). For example Lavoie (1974) applies a filter to the topography to prevent the propagation of spurious waves within the model. Another more desirable solution is to use a parameterization of these sub-grid scale fluxes that extracts energy from the averaged equations in a manner that is consistent with reality (Pielke,1984). Thus it can be seen that numerical methods used to solve the governing equations can directly affect the choice of model formulation. It has been shown that numerical models are limited in their ability to model the real atmosphere by their inability to resolve the vast spectrum of real length and time scales. In addition, limitations are also invoked by the extent to which the physical processes in the real atmosphere are represented by the model formulation. The increase in resolution or the inclusion of more physics in the model by the use of more general formulations of the governing equations, leads to a large increase in computing power required. This level of computing power may not be available for technological or financial reasons and as such has tended to limit models toward being more phenomenon specific. 1.3 Choosing a Model At this point, with a better idea of some of the limitations involved, a look at the phenomena of interest is in order. This needs to be done within the context of striking a - 1 1 -balance between including the physics required to simulate the phenomena of interest and keeping the model simple and economical. As stated earlier, a one layer model appears to be a good choice for a starting point. The model must solve the primitive equations which include the thermodynamic and momentum equations near the surface. Because the model is one-level and does not calculate vertical velocities, it cannot be mass consistent. For use in this study the model must have or present the possibility for the inclusion of spatially variable surface roughness and heat exchange between the earth's surface and the atmosphere. The model is to be used over periods of time long enough for the surface energy balance to change substantially due to diurnal variations in solar radiation inputs. Therefore the model must also accommodate the inclusion of temporally varying thermal forcing in order to simulate the diurnal cycle of heating and cooling. The model chosen for this study was one originally formulated by Danard (1977). The Danard model was subsequently modified by Mass (1981) and Mass and Dempsey (1985a). This version of the model is described in detail by Mass and Dempsey (1985b). A copy of this model was very kindly made available by Mass and Dempsey for use by the author in this study. The Mass and Dempsey (1985b) model distinguishes between land and sea surface characteristics, but does not resolve spatial variations over the land surface. The model was subsequently modified by the author in order to enhance its capability to deal with spatially varying surface roughness and thermal characteristics as well as temporally varying thermal inputs. The model was also modified in such a way as to make the surface friction terms in the equations of motion dependent upon stability. These modifications are described later in the text. For the sake of completeness there are a number of details presented in this thesis that are also presented in Mass and Dempsey (1985b). There are also a number of model details that are not described within this thesis but do appear in Mass and Dempsey (1985b). Where not specifically stated, details of the - 1 2 -model should be assumed to be unchanged from the version of the model presented by M a s s and Dempsey (1985b). The work of the author in this study should be considered an extension of work performed by M a s s and Dempsey in the hope that an already useful tool might be made even more useful by a s m a l l number of modifications. -13-2. THE MODEL The following chapter is intended to explain the workings of the model used in this study, and some of the more general aspects that are common between models of a similar nature. Throughout this chapter there are references to appendices made with the intent of maintaining an acceptable level of continuity in what can be a very complicated explanation. Assumptions used in the formulation of the model are stated, and with the exception of those cases where it lends to the explanation, the implications of these assumptions are left to be examined elsewhere in the thesis. The general structure of the model is in groups and sub-groups of commands which act upon variables, therefore the clearest way to structure the explanation is in a pattern which follows the program structure. Before this can be done, however, there must be an explanation of the general model formulation and the physical principles underlying the model equations. 2.1 Model Formulation The model makes use of a piecewise linear temperature structure to parameterize the vertical temperature profile throughout the model domain, (see Figure 2.1.1) The temperature (Tr) and height (Zr) of a reference pressure level are assumed known. The temperature profile is then extrapolated down to the top of the layer of topographic influence (Z^), using an input lapse rate which is calculated from real atmospheric conditions. The top of the layer of topographic influence follows the ground surface at an elevation which is a distance H above the ground surface. Below this level, the temperature structure is a linear interpolation between the temperature at the top of the layer of topographic influence (TQ) and the temperature at the surface (Tg). The surface temperature is allowed to vary due to diabatic and adiabatic effects. Changes in temperature have the effect of changing the interpolated temperature structure between V E R T I C R L T E M P E R R T U R E " S T R U C T U R E P R R R M E T E R I Z R T I O N • r-n F i g u r e 2 . 1 . 1 V e r t i c a l T e m p e r a t u r e S t r u c t u r e P a r a m e t e r i z a t i o n -15-the surface (Z ) and the top of the layer of topographic influence. The temperature structure of the layer between Z n and the reference level is held constant through time. This temperature structure over the height Z f to Z g is vertically integrated to give pressure at the surface. This results in a surface pressure that is dependent upon the vertical temperature structure which, in turn is only dependent upon the modeled temperature at the surface. It is assumed here that the mechanical and thermal effects of the surface are present only up to a certain height in the atmosphere. The counterpart to this layer of topographic influence in the real atmosphere is the planetary boundary layer. It is recognized that the level of thermal effect and the level to which the atmosphere is affected by mechanical forcing are not the same. Further to this, these levels are known to change through time and space. These spatial and temporal variations in the boundary layer height are not modeled explicitly and as such the model is subject to errors when the vertical structure of the model diverges from the physical reality of the atmosphere. The above parameterization uses the pressure (Pr) and the temperature as boundary conditions at the top of the domain. This is done by choosing the upper boundary of the domain to be a surface of constant pressure and specifying the elevation and temperature over that surface. For this study, the reference pressure surface used is the 850mb surface. These upper boundary conditions specify the synoptic scale forcing over the model domain. It is clear at this point that through the vertical integration, the pressure at any one point on the ground surface is dependent upon the temperature and therefore the density structure overlying that point. If the temperature of the layer above the uppermost extent of topographic influence, Z n to Z r , is constant through time, then using this parameterization, variations of surface pressure with time at any point must be the result of advective changes in temperature, fluxes of heat to or away from the surface itself and diffusion of heat from surrounding areas. - 1 6 -The model is set in the form of two second order partial differential vector equations that describe the tendency of temperature and wind velocity at the surface. ! l i = _v s.v o T s + Ei{!l^s + v s - V a l n P s } + £ _ + K T V H 2 T S o t C p o L Cp ! l l = - V S - V 0 V S - f k x V s - { g V a Z s + R T s V 0 l n P s } + F + K M v a 2\7 s 3 t In the above equations, T g is the surface temperature, P s is the surface pressure, Q is the input heat energy, K t and Km are the coefficients of diffusion for heat and momentum respectively, R is the universal gas constant (2S7 J K ' ^ kg"1) , c is the specific heat of air at constant pressure (1004 J K"1 kg" 1), V"s is the surface wind velocity vector, Z g is the surface elevation, f is the Coriolis parameter (1.10x10"^ s" 1), F is the force due to friction and g is the acceleration due to gravity (9.8 m s" ). The above equations are in sigma coordinates. The transformation to sigma coordinates involves the transformation of the vertical coordinate, z in this case, to a log-pressure coordinate. The value of sigma at a point at or above a point on the surface is the logarithm of the pressure at that point divided by the logarithm of the pressure at the surface underlying that point. This means that sigma has a value of unit}' everywhere on the surface. This is of particular significance in this case because the model is one-level and that level is the surface. In effect, sigma coordinates are terrain-following coordinates. Appendix I shows the effect of the sigma coordinate transformation on the vertical coordinate and points out how it is used in the formulation of equations 2.1.1 and 2.1.2. The full derivation of these equations is presented in Appendix II. The first term on the right-hand side of the temperature tendency equation (eq.2.1.1) is a term that accounts for the advection of temperature. The second term on the right-hand eq.2.1.1 eq.2.1.2 -17-side is one that represents the change of temperature due to a change in pressure, for example, adiabatic warming or cooling. The next term is the diabatic forcing term and represents the flux of heat to or from the surface at any one point. The last term accounts for the diffusion of heat to or away from a point. This term also helps to suppress numerical instability in the scheme used for solving these equations (Pielke,1984; Anderson, Tannehill and Pletcher,1984; Haltiner and Williams, 1980). Equation 2.1.2 is the velocity tendency equation. The first term on the right-hand side is a momentum advection term. The second term is one that represents the effect of a rotating reference frame on the motion of the air within the domain. The model uses a constant value of f. It is assumed that the Coriolis parameter (f), which in reality varies with latitude, is constant within the domain of the model. This is justified because the domain is of little latitudinal extent, thus f changes by less than 2% over the north-south dimension of the domain. The next term in the equation is the pressure gradient term. This term will eventually be manipulated to be a function of the surface temperature through the vertical integration of the temperature structure. The fourth term is a term that represents the force of friction on the air as it moves over the surface. The final term in the equation accounts for diffusion of momentum. This has a similar role to the diffusion term in the temperature tendency equation in that part of its purpose is to damp out numerical instabilities in the solution of the model equations. As is the case in the temperature tendency equation, the diffusion term acts to damp out instabilities by reducing large gradients created by spurious high wavenumber perturbations. The manipulation of the grad-Pg term used in both the temperature and velocity tendency equations is presented by Mass and Dempsey (1985b). In this manipulation any explicit reference to pressure or its gradient at the surface is removed. This transformation casts the pressure at the surface in terms of the vertical temperature structure which, through the parameterization used, is a function of T r , T n and T g . Once this manipulation - 1 8 -is carried out, the tendency equations are coded in a finite difference form to be solved using a computer. The model uses an Arakawa "C" staggered grid as described by Mesinger and Arakawa (1976). The advection terms in the momentum equations are centre differenced in space (Gerrity et aL,1972) and the momentum diffusion terms are differenced using a scheme by Danard (1971). The centred differencing of the temperature advection terms and the temperature diffusion terms are according to Danard (1971). 2.2 Model Equation Solution 2.2.1 Initialization The tendency equations used in this model need some initial conditions. The fields of temperature and wind velocity must have an initial starting point from which the integration can progress through time. Once the initialization has taken place, the model can be forced diabatically by fluxes of heat or by changes of friction which, through a stability correction, are a function of heating. There are two ways of initializing the wind and temperature fields. The first method is a first guess and subsequent integration to steady state. The model extrapolates the upper lapse rate down to the surface and from that temperature structure calculates the surface pressure and subsequently the surface pressure gradient..This pressure gradient is used in the balance between friction, the pressure gradient and the Coriolis "force" to obtain a first guess for the surface windfield. Together with the surface temperature field which has already been calculated at this point, the model has a complete set of first guess fields. From this point the model must be run to steady state with no diabatic heating. This has the effect of bringing all of the fields to a state where they are consistent with the physics of the model. The term steady state as it is used here refers to a condition in which none of the fields are changing with time. Inherent in this method for solving the equations of motion for the - 1 9 -atmosphere, is a certain amount of numerical error due to the limits of precision of the computer. This means that the fields will never reach a state of exact stationarity. The method used in this case, is to calculate the average tendency of the velocity field, and when that average tendency reaches a sufficiently low level the field is assumed to be stationary. This level of stationarity is a variable that is declared in the input of the model run and is usually 10"^ m s"^ . The second method of initializing the model is to "hotstart" the model using previously modeled fields as a starting point for the integration. This method is used in order to eliminate the cost of integrating to steady state each time the model is run. As an example, if one wanted to make three different runs with varying model parameters, the three runs could be made from the same set of initial fields without the cost in time and computer resources of three integrations to steady state. It is important to note here that, because the fields must be consistent with the physics of the model, there must be no step changes in any of the controlling parameters such as diabatic heating or friction. For example, factors such as diabatic forcing must be kept constant from before to after the "hotstart" procedure. If one or more of these factors is changed more than a small amount over the "hotstart", the model must respond to what amounts to a step discontinuity. Because the model uses finite differences to represent the governing equations, a step change cannot be properly handled. This inability to handle step changes can cause spurious results to appear in the solution or cause the failure of the model run. 2.2.2 Evolution through time of modeled fields Once the model has a starting point made up of consistent surface fields of wind velocity and temperature, the model must then use the tendency equations to determine the state of the fields at some time in the future. This process is called time stepping. The temperature tendency equation calculates the change in temperature with time over the domain due to diabatic and adiabatic effects. Using the temperature tendency -20-equation and the present wind and temperature fields, the change with time of the surface temperature can be calculated. Using this rate of change with time, the value of the surface temperature can be extrapolated to the next time step. On the first step the extrapolation is made linearly, with subsequent steps being made using a modified second order Adams Bashforth method (Mass and Dempsey, 1985b; Burden, Faires and Reynolds, 1981). Once the updated temperature field is calculated, it is used in the velocity tendency equation to calculate the derivative with respect to time of the velocity over the domain. The velocity field is then stepped forward in time in a similar way to the temperature field, to arrive at an updated velocity field. Thus, one time step involves the calculation of a new temperature field and a new velocity field. Strictly speaking, because the above method uses the already updated temperature field to calculate the velocity tendency, the equations are not solved simultaneously. This is of little consequence since the tendency equations do not involve terms that change rapidly with time or have within them severe non-linearities. The length of the time step used in the time stepping process is of great importance to the accuracy of the fields being modeled. If the time step is too great this can lead to numerical instability. This numerical instability is the result of the inability of the finite difference scheme to represent exactly, pertubations in the modeled fields. This inability creates numerical dispersion and diffusion specific to the equations being solved and the finite difference schemes used. Numerical instabilities can cause perturbations to grow undamped and overwhelm the desired solution. Some of the effects of truncation resulting from the use of finite difference representations are described by Gerrity et al. (1972). Pielke(1984), Anderson, Tannehill and Pletcher (1984), Haltiner and Williams (1980) and Holton (1979) all provide treatments of numerical instability in various differencing schemes. -21-Using a time step that is too short, although keeping the solution stable, ma}' result in less than optimum use of computer resources and poor performance by the model. For a given time frame to be modeled, it will take a number of time steps equal to the length of time modeled divided by the length of the time step. It is clear that the shorter the time step is, the more steps that will be necessary to model a given period of real time. This results in the use of relatively greater amounts of computer resources. Aside from the uneconomical use of computer resources, there is another drawback in the use of a time step that is too short. The method for calculating any of the new fields involves a great number of calculations. Each one of these calculations is subject to roundoff error. This roundoff error is the result of the floating point characteristics of the computer. In effect the actual numbers being used in a floating point operation must be rounded off to as many decimal places as the floating point representation of the computer will allow. In the case of the Floating Point Systems AP164/MAX array processor which was used for most of the modeling in this study, the floating point numbers are represented to sixteen decimal places. The error induced by this sort of truncation may, upon first inspection, seem very small but when one considers the number of floating operations needed to go through one time step, the roundoff error becomes more significant. It is for this reason that the time step should be kept as long as possible without adversely affecting model stability and, therefore, the modeled solution. 2.3 Model Outputs The input/output routines of the model have been written in such a manner as to facilitate the output of any or all of the model arrays. The arrays can be output at any interval which is an integral number of time steps in length. The values in these arrays are the instantaneous values of the model variables at the time of output and are in no way temporally averaged. These arrays can be output during the integration to steady - 2 2 -state or during subsequent integration using diabatic forcing. This output facility was used a great deal in the development stages of the modeling. The validation of a model of this nature will in most cases be done with temporally averaged windfield data. For this reason, the model has written into it the ability to output hourly averaged windfields. The averaging operator is a simple mean of the fields of U and V components of the modeled windfield. The averaging is done at the end of each hour and the averaged values are output to a file which is separate from the above mentioned output arraj's. -23-3. SURFACE PARAMETERIZATION By its nature, the method of finite-differences, when applied to the solution of differential equations, carries with it the discretization of what are in reality continuous fields. This discretization defines spatial and temporal resolutions which in turn determine minimum length and time scales of features which can be represented by the solutions. The atmosphere, unlike its discretized analog, is continuous and obeys physical principles which dictate spatial and temporal scales well beyond the limitations imposed by the finite-difference representation of scalar and vector fields. As an example of this "truncation of scale", the resolution of the model used in this study was approximately four kilometres in both of the spatial dimensions and sixty seconds in the time dimension. Clearly there are many atmospheric phenomena which operate at length and time scales much smaller than this. Modeling of atmospheric phenomena is, to a great extent, the modeling of the flow of energy from one point to another as well as the cascade of energy from one scale to another. This cascade, with very few exceptions, is from the larger scales to the smaller scales in the real atmosphere. How then, does a model having lower limits of space and time scales that are much larger than the lower bounds of the real atmosphere, represent this cascade? This question seems particularly pointed when one considers the fact that there are at least five orders of magnitude in space and three orders of magnitude in time between the model resolution mentioned above and small scale turbulence in the real atmosphere (Atkinson, 1984; Oke.1978). (see Figure 3.1) The answer to this question lies within surface parameterization or parameterization of sub-grid scale fluxes. As the term implies, the parameterization of sub-grid scale fluxes is the approximation of this cascade of energy from larger scales to smaller scales in both length and time. This model makes use of two main parameterizations. The parameterization of surface friction which is the approximation of the downward flux of momentum near the surface and the 1 C M 10 C M 1 M 10 M 100 M 1 K M 10 K M 1 0 2 K M 1 0 3 K M 1Q 4 KM 1 0 5 K M M O QUASI-YR BIENNIAL OSC'NS INDEX CYCLES WK D A Y HORIZONTAL INERTIAL OSC'NS HR VERTICAL BUOYANCY MIN OSC'NS 10 10" 10" 10" 10" Charac te r i s t i c hor i zonta l s c a l e L ( m ) i I Figure 3.1 S p a t i a l and temporal c h a r a c t e r i s i c s of atmospheric phenomena ( a f t e r Smagorinsky 1981) . -25-parameterization of diabatic heating/cooling which is the estimation of the flux of sensible heat to or from the atmosphere within the boundary layer. 3.1 D i a b a t i c H e a t i n g / C o o l i n g P a r a m e t e r i z a t i o n Mass and Dempsey (1985b) use values for the diabatic forcing which are specified by heating/cooling curves. These are designed to be representative of the rates of heating and cooling throughout one diurnal cj^cle. These curves differ between land and ocean surfaces but not within those two tj^pes of surfaces. It is the intention of the present study to modif}' the diabatic forcing parameterization so as to make it more realistic as well as to make it handle spatial^ varjang surface heating and cooling. With this in mind, a new parameterization was devised. One of the drawbacks of the one-level model is that, where multi-level models can use vertical temperature gradients to parameterize transport of sensible heat from the earth surface to the atmosphere, the one-level model cannot. Therefore there can be no closure of the surface energy budget. This means that surface energy budget terms must be found elsewhere. The most logical choice lies in the parameterization of surface heating through radiative inputs. The parameterization used in this study is based on the net all-wave radiation balance derived from real data for the shortwave term and a simple model for the longwave term. Using measured values of direct and global solar radiation, the value of the diffuse component of the solar radiation was calculated. These values were subsequently modified to account for the spatially varying slope angle and slope azimuth over the domain. The correction to the direct beam radiation is made using a formulation presented by Hay and McKay (1985) for the irradiance on a sloping surface: S = I cos(i) eq.3.1.1 - 2 6 -where S is the direct irradiance for a sloping surface, I is the radiant intensity at normal incidence and i is the angle of incidence between the Sun and a normal to the surface. A similar correction to the diffuse component for slope angle and azimuth is made using an anisotropic slope irradiance model presented by Hay and McKay (1985): D = 0{KCOS(i)/cos ( z ) + 0.5(1.0-<)(1.0+COS(a)} eq.3.1.2 where D g is diffuse solar irradiance for an inclined surface, D is the diffuse sky radiation, < is the anisotropy index, z is the zenith angle and a is the slope angle. Details of the calculations involved in the correction for slope and zenith angle to the direct and diffuse beam solar radiation appear in Appendix IV. The sum of the direct and diffuse components is then multiplied by 1.0 minus the surface albedo to obtain the net shortwave component of the budget. The value used over the whole domain for the albedo was 0.2 (Steyn,1980). The longwave component of the budget was handled by the Idso and Jackson (1969) formulation for the net longwave radiation which makes use of the near surface air temperature: L* = aT a4{-0.216exp[-7.77xl0-4(273-T a)2]} eq.3.1.3 where L is net longwave radiation, a is the Stefan-Boltzmann constant (5.67x10 W m K"4) and T a is the near surface air temperature. The short and longwave components were then added together to arrive at a figure for net radiation. - 2 7 -The sensible heat term was calculated using a Bowen ratio method. The formulation of this is as follows: 0 h = f3(Q* - Q s )/((1+B) eq.3.1.4 where $ is the Bowen ratio (the ratio of the sensible heat to the latent heat fluxes), Q* is the net all-wave radiation and Q s is the storage term. The storage fraction used for the domain was 0.2 from Oke et al. (1981). This value was held constant through time and was the same over the entire domain. Although this value is known to change with time and surface characteristics (Oke et al., 1981), in the absence of more data concerning the surface characteristics, one value thought to be representative of a regional average was used. The values of the Bowen ratio used over the domain did not vary in time but were significantly different from point to point within the domain. Cleugh and Oke (1986), in a study involving an area within the model domain give typical values of .46 and 1.28 for the Bowen ratio in rural and suburban areas respectively. Values of the Bowen ratio based on the relative proportions of suburban and rural land use within a grid square were used, (see Figure 3.1.1) This thesis makes use of two modeled days during the summer months of 1985. Figures 3.1.2 and 3.1.3 show the modeled surface energj' budget terms through these two days. The budget figures presented are representative of the modeled surface energy budget at the grid node closest to. a research facility operated by Dr. T.R.Oke of the Department of Geography at the University of British Columbia. This site is in a suburban area of South Vancouver and is described by Cleugh and Oke (1986). Figure 3.1.4 shows tj'pical measured values for budget terms taken at this site. The values plotted in Figure 3.1.4 are averages taken over the months of June through September 1985. As can be F i g u r e 3 . 1 . 2 Modeled E n e r g y Budget T e r m s , J u l y 20,1985 F i g u r e 3 . 1 . 3 Modeled E n e r g y Budget T e r m s , A u g u s t 23 ,1985 Figure 3.1.4 Measured Energy Budget Terms (After Cleugh and Oke, 1986) CHRNGE IN Ps vs HEIGHT OF TOPOGRRPHIC I N F L U E N C E F i g u r e 3.1.5 S e n s i t i v i t y of S u r f a c e P r e s s u r e Change to H e i g h t of Topographic I n f l u e n c e . T h i s p l o t shows the change of p r e s s u r e at the s u r f a c e i n response to an i n c r e a s e i n temperature at the s u r f a c e . The temperature change at the s u r f a c e ranges from IK to 10K with the d i a g o n a l l i n e s r e p r e s e n t i n g the c o r r e s p o n d i n g change i n pressure at the s u r f a c e f o r g i v e n h e i g h t s of t o p o g r a p h i c i n f l u e n c e . - _ 3 3 -seen, the instantaneous modeled values are slightly higher than the averaged values but the pattern of diurnal variat ion is very well modeled. In summariz ing Budyko's results, Oke (1978) gives a representative value for over water of 0.1. This value is multiplied by a factor of 0.1 in order to account for the near constant sea surface temperature dur ing the diurnal cycle and the resultant limited vertical mix ing of sensible heat (Oke,1978). A s can be seen in Figure 3.1.5, the sensitivity of surface pressure to temperature variations is l inearly related to the depth to which the heat from the surface is allowed to m i x . Because the mixed layer height over the sea surface is smal l , the response of the surface pressure to solar inputs at the surface is expected to be smal l as wel l . Once a value of Q n is calculated, the heat is assumed to be vertical ly mixed through the layer of topographic influence. The m i x i n g is assumed to spread the heat in a l inearly decreasing (upward) fashion, up to the height of topographic influence, above which no heat is added. The value of is used to calculate the surface temperature tendency. The derivation of this formulation is presented in Appendix V. The formulation is consistent wi th the modeled l inearly decreasing (upward) distribution of sensible heat through the layer of topographic influence. 3.2 Surface Friction Parameterization The atmosphere can be thought of as having three mechanisms for transporting momentum. These three mechanisms act at different scales ranging in length from the global scale to the molecular scale. This tremendous var iat ion in length scales is = RTsQh P sc pH eq.3.1.5 Equation 3.1.5 gives the formulation for the surface temperature tendency given Q ^ . -34-accompanied by an equally large spread of time scales through which these mechanisms act. The first of these mechanisms is the horizontal or vertical propagation of coherent waves in the atmosphere. These waves can range in size from metres to thousands of kilometres as in the case of planetary waves (Beer, 1974). In terms of friction parameterization, the only waves of significance are the vertically propagating waves. Because the model used is one level, it cannot represent vertically propagating waves. This has importance not only in terms of the pressure perturbation at the surface caused by these waves but in terms of the energy transported to or from the surface. The inability to model these waves makes the use of the model inappropriate where waves of this nature would make a significant contribution to the vertical fluxes at the surface or the variation in the surface pressure field. The remaining two methods of transporting momentum in the atmosphere are turbulent motion and molecular diffusion. The scale of these motions ranges from tens of metres to the molecular scale (Atkinson, 1984). These motions are stochastic in nature and are normally quantified using statistical methods. The transport contribution of molecular diffusion is very small when compared to that of turbulent motions in the atmosphere (Pasquill and Smith, 1983). Therefore, when subsequent reference is made to turbulent transport, contributions from molecular diffusion are assumed to be included. In order to quantify the force of surface friction upon a unit volume of atmosphere, a link between friction and vertical turbulent transport of momentum must first be established. This link is most easily seen by examining a unit volume of air near the surface which has no net vertical motion and no net vertical forces acting upon it. Horizontal momentum is being transported by turbulent eddies into the unit volume from the top and out of the volume at the bottom. Transport of horizontal momentum vertically across a horizontal plane is equivalent to a horizontal stress over that plane. The difference between the stress over the two horizontal surfaces bounding the layer is the -35-divergence of shear stress. This divergence of stress is equivalent to a force per unit volume acting in the opposite direction to the mean flow. This frictional force needs to be parameterized in order to make it usable in the equations of motion employed in the model. The method used by Mass and Dempsey (1985b) is one in which the stress is parameterized by: S S = pC r jV s |V s | - eq.3.2.1 where S g is the horizontal shear stress near the surface, is the coefficient of drag and V"s is the wind velocity vector. If a linear profile of stress is assumed, where stress vanishes at the top of the layer of topographic influence, then: P = - l f S h - S s i 3 - C d V S l V S | where F is the force per unit volume exerted upon that volume, is the fluid density, S g and S n are the shear stresses at the top and bottom of the layer respectively and H is the depth of the layer. As can be seen in equation 3.2.2, the friction force on a unit volume of air is inversely proportional to the thickness of the boundary layer. The modeled value of H is constant through time, therefore the friction is multiplied by a factor to account for the variable height of the frictional boundary layer. Mass and Dempsey (19S5b) use values of 2.0 and 4.0 for Aj in the day and night respectively. These figures are specific to the model runs * used in that study. In the present study the value of A- is linked to the value Q , the net radiation at the surface. Steyn and Oke (1982) present values for the inversion height (Zj) through a summer day over suburban Vancouver. The thickness of the inversion layer ranges from as low as -36-50m at night to a maximum of 400m near mid-day. The modeled value of H is 1000m. therefore the friction must be multiplied bj' a factor which reflects the ratio of the modeled to the real inversion heights. This diurnal pattern is accounted for in a very simple fashion by specification of Aj as shown in equation 3.2.4. This adjustment acts over the whole domain. Mass and Dempsey (1985b) increase the stress by a factor of 2.8 under neutral to stable conditions. This is suggested by Deardorff (1972) to account for the fact that the stress vanishes at a height lower than the inversion height under these conditions. The same procedure is followed in this study. This gives: F = -Aj2.8C dV s[V s| e q. 3 - 2 3 H where: Ai = 3.0 - ( Q * / 1 0 0 ) Q* < 100 Win"2 = 2.0 0* > 100 Wm-2 eq.3.2.4 Mass and Dempsey use the above formulation with coefficient of drag (C^) specified as one value over land and one value over water. The above scheme can be used with a spatially varying, constant in time over the domain in order to accomodate spatially varying surface roughness.. Instead, equation 3.2.3 was used in conjuction with a method of specifying that can handle not only spatially varying surface roughness but also make stability dependent. This method uses the formulation for used b}r Keyser and Anthes (1977): k 2 C d =, 7-z eq.3.2.5 ( l n ( z / z 0 ) - ^ m ( z / L ) } 2 where: » m = i . : i n [ r i - 1 6 ( z / L l l 2 / 3 * [ 1 - 1 6 ' z / L ) ] 1 / 2 * x \ - 1 . 7 3 ( t a n - M 2 [ ' - 1 6 ( ^ ] 1 / 3 * *) - L O S ) eq .3 .2 .6 L = - u * 3 p C P T eq.3.2.7 kqQh In equations 3.2.5 to 3.2.7, k is the von Karman constant (.40), z n is the roughness length, L is the Monin-Obukhov stability length, u* is the friction velocity, is the sensible heat flux, z is the height above the surface, T is the temperature at height z and ij> m is a dimensionless function of stability. Spatial variabilit}' of surface roughness is added to the model by the use of a spatially varying value of ZQ. In order to specify spatially varying surface roughness, values of ZQ taken from literature were assigned over the domain, based on estimates of surface roughness from air photographs of the domain (Pielke,1984; Stejm,1980; Silversides, 1978; Oke, 1978; Counihan, 1975). A map of the roughness length over the domain is presented in Figure 3.2.1. In the above formulation, a stability dependence is brought into the formulation through the use of a logarithmic wind profile that includes the non-dimensional stability function. In a modified formulation, using: C d = JvTp eq.3.2.8 equation 3.2.5 can be expressed as: * 2 a k 2 | V s | 2 . [ln(z/z0) - *m(z/D]2 eq'3'2-9 CO 00 Figure 3.2 .1 Roughness Length Over the Modeled Domain -39-This leads to the non-linear equation: U* = eq.3.2.10 [ l n ( z / z 0 ) - 7m{. •}] Keyser and Anthes (1977), in a somewhat more complex model solve iteratively for u* possibly large number of iterations at each time step at each grid node. This would involve a great number of extra calculations at each time step to include stability dependence in the model. This also assumes that the method used (fixed point iteration, Newton's method, bisection, secant, etc.) will converge quickly from some starting point. In this assumption lie two problems; the first is that the best one could hope for with a well conditioned problem (which this is not) is second order convergence. This would still require a very large number of calculations at each time step. Secondly, the starting point for the iteration cannot always be chosen to give the fastest convergence or to even guarantee convergence at all. (Burden, Faires and Reynolds, 1981) In order to get around this problem, an approximation was made. Instead of using an iterative technique within the model itself, the iteration process was brought outside the model by the use of a statistically defined function for u*. Values of mean wind speed (U), roughness length (ZQ), temperature (Tg^ and the sensible heat flux (Q n) were put into equation 3.2.10 and an iterative solution for the stability corrected u* was found. This was done using different values of the above variables in permutations chosen to include a wide range of stabilities (-.001>(Z/L)>-15.0). The ratios of the stability corrected u* to the neutral u* were calculated for all cases and a multiple linear regression was performed in order to find a function for the correction factor. The and C j through equation 3.2.8. The iterative solution of this equation would involve a -40-regression was done with various combinations and transformations of the independent variables in equations 3.2.5 to 3.2.7. Several trials were needed to find a parameterization that was not only reasonably accurate but smooth near neutrality. The best results were obtained when the independent variables chosen for use in the regression had values of zero at neutrality or varied slowly about zero at neutrality. This was achieved by choosing groups of independent variables that, where possible, were multiples of the surface heat flux. This seems reasonable in light of the fact that with no net heat flux at the surface, the stabilit}' of the atmosphere approaches neutrality. Numerically, this can also be seen as the Monin-Obukhov stability parameter goes to infinity as Q n goes to zero, forcing the value of the stability function used in equation 3.2.4 to zero. This has the effect of smoothly forcing the statistically defined correction function to a value near 1.0 at neutrality. It was found that the correction for stability in the stable case was small compared to the correction in the unstable case. In order for the. atmosphere to sustain stable stratification, the wind speeds and therefore u* must be very low. The friction 9 parameterization is a function of u* , thus the effects of any correction to u* at very low wind speeds is quite small. Therefore, in the interest of simplicity, only one function (eq.3.2.11) was used. This function does adjust u* a small amount in the proper direction for the stable case but the accuracy of the adjustment is highly questionable. It was thought that the small loss in accuracy of the parameterization was worth the gains involved with making the scheme more robust numerically by making the function smooth near neutrality. The function used to correct u* for stability was: u(unst) = u * ( n t r l ) [ 1 . 1 6 9 + .043741 n(-z/L) + .002517{ -^U] eq.3.2.ll - 4 1 -Figure 3.2.2 shows the performance of equation 3.2.11 in producing correction factors in lieu of the more computationally intensive iterative process. The ordinate axis is the value of the u* adjustment arrived at by the iterative (exact) process and the vertical axis is the value of the adjustment arrived at by the regression (approximate) formula. The statistics of the regression are also presented in Figure 3.2.2. The values of A r and B r are, respectively, the slope and the intercept of the ordinary least squares regression line between the values arrived at using the iterative process and those arrived at using equation 3.2.11, at various stabilities. SQ and Sp are the standard deviations of the exact and approximated u* correction factors. Also presented are the root mean square error (RMSE), the systematic root mean square error (RMSEg), the unsystematic root mean square error (RMSEU) and the index of agreement. The value of the index of agreement ranges from 0.0 in the case of no agreement to 1.0 in the case of perfect agreement. The formulations used to calculate these indices are presented in Appendix VI. The value of the u* stability correction factor calculated using equation 3.2.11 through the modeled days is presented in Figures 3.2.3 and 3.2.4. These values are output at the same grid point as the energy budget terms in the previous section. Though there seems to be no way to quantitatively assess the validity of this adjustment of u* for stability, the curves appear to be reasonable. In the early morning and late evening when sensible heat flux is toward the ground, the atmosphere is stable and the turbulent transport of horizontal momentum is suppressed (Oke, 1978). During this period the adjustment acts to reduce u*. During hours with high solar radiation loading, large fluxes of sensible heat are transported into the atmosphere, causing instability near the surface. This instability tends to increase the turbulent intensity near the surface and thereby enhance the downward transport of horizontal momentum. This increase in shear stress appears in the u* stability correction as a peak near mid-day. Figure 3.2.2 U* S t a b i l i t y Correction - Exact vs Approximated -e<7-- 7 7 " 4. EXPLORATORY MODELING The model being presented embodies a number of changes from the original model obtained from Mass and Dempsey (1985b). The changes made to the model are for the most part in the parameterization used. The original model resolved differences in surface roughness and thermal characteristics between land and sea surfaces but not within those categories. The original model also used a somewhat simplified parameterization for the diurnal pattern of heating and cooling over both land and sea surfaces. The model as used for runs presented here, uses a parameterization scheme that accounts for spatial variability in surface roughness and heating. This modified formulation also links the temporal pattern of heating and cooling to a parameterized surface energy budget. In order to evaluate the effects of these changes, a modeling strategy that examines the behavior of the model in response to the new parameterization was necessary. 4.1 Testing of Surface Friction Parameterization At the upper boundary of the model domain, the effects of friction are assumed to be negligible. At this elevation the steady state flow is geostrophic and runs parallel to the contours of the 850mb surface. Below this level, in both the real and modeled atmospheres, the balance of forces changes to include the force of friction. This frictional force is the result of shear stress and has the effect of creating a cross-isobaric component to the flow. Over homogeneous roughness, the flow is parallel to the geostrophic flow at the top of the frictional boundary layer and in the Northern Hemisphere, backs with vertical distance below that level. At a level just above the surface, the flow has a reduced speed and direction which are determined by the roughness qualities of the surface. This rotation of the flow with height is called the Ekman spiral, named for the Swedish oceanographer V.W. Ekman who first -46 derived a solution to the analogous situation in the surface layer of the ocean (Holton, 1979). This analytic solution is presented in many places in the literature in many forms (Haltiner and Williams, 1980; Holton,1979; Pielke,1984) These analytic treatments of the subject are helpful in understanding the general form of the Ekman spiral and its relation to the height of the frictional boundary layer. These forms of the solution involve the use of a parameter which has within it a value for the vertical diffusion coefficient. In the one layer model, only the horizontal coefficient, of diffusion is specified. Therefore the analytic solution can only be of use in a qualitative way to assess model response to surface friction parameterization. In assessing the performance of the friction parameterization, it is clear that modeled results must be compared in a quantitative way with well understood behavior of the frictionally affected layer of the atmosphere. The parameters chosen for comparison must be common between modeled results and values calculable using current knowledge of the atmospheric boundary layer. Assessment of the model's friction parameterization was done by making model runs over fiat terrain with homogeneous roughness. This was done in order to determine the modeled value of the angle between the surface flow ( B ) and the geostrophic flow at the top of the boundary layer and the magnitude of the wind speed just above the surface (U*/V ). Given U*/V and 8 , the accuracy of the modeled balance between the Coriolis force, the pressure gradient force and the frictional force can be determined. Estoque (1973) presents in a graphical form, values of U*/V and 6 as functions of surface Rossby number (V /fZg). Panofsky and Dutton '(1984) give a formulation for the angle as:. tan(B) = -l n ( Z 1 / Z 0 ) - A eq.4.1.l where Z- is the thickness of the boundary laj'er. -47-Using this formulation and values suggested by Panofsky and Dutton of 0.0 and 5.0 for A and B respectively, B can be calculated given the surface roughness length and the height of the boundary layer. The values for A and B seem to carry with them a certain amount of arbitrariness but agree quite closely with the Estoque (1973) values in sample calculations. It is also assumed in the calculation that Z>100Zn. Panofsky and Dutton also present in a modified version of an unpublished figure of H. Lettau, a nomogram for U* as a function of V and 1/Znf. This allows the calculation of U : | :/V . Several model runs were made to obtain values of B and U*/V . Each run was made with a different roughness length ranging from .00021 m to 2.00 m. The geostrophic wind at 850mb (approximately 1500 m) was from 270 degrees at 10 m s"*. Table 4.1.1 gives modeled and literature values for B and U : !7Vg at various roughness lengths. As can be seen values of B and U*/V given by Estoque agree quite closely with values given by Panofsky and Dutton. The modeled values of B and U*/V show good agreement with these two sets of values up to about ZQ=.5m.. At roughness values greater than this, the model seems to overestimate the turning and the magnitude at the surface by as much as 30 percent, (see Figure 4.1.1) This result is not unexpected given that the modeled elevation is 10m, giving the smaller Z/ZQ values very much lower than those used in the assumption that Z/ZQ>100. Given a modeled height of 10m, the comparison is only valid in the strictest sense for roughness lengths of .10m or less. The only test that falls within this limitation is that with a roughness length of .00021m. The modeled results, although divergent from theoretical values at high roughness, follow the theory very closely at roughness values closer to the .10m limit. The model was also run over a step discontinuity in roughness to observe the distance needed for the flow to become fully adjusted to the change in friction. This was done so that the flow was at right angles to the discontinuity. In the case of a very large jump in roughness, it was found that the flow was not fully adjusted until it had reached a point at least five grid spacings downstream from the discontinuity. This implies that in the case of -48-FRICTION PARAMETERIZATION COMPARISON Modeled Estoque Panofsky and Dutton z 0 u * / v g u * / v g u * / v g .002 .027 .027 .030 .25 .043 .042 .040 .50 .049 .044 .044 1.06 .057 .046 .046 1.25 .058 .047 .047 2.00 .062 .048 .050 Z 0 .002 .25 .50 1.06 1.25 2.00 Beta 15.5 27.5 29.6 43.7 44.3 53.5 Beta 16.7 30.7 32.2 35.0 35.2 36.5 Beta 18.0 31.1 33.3 36.1 36.8 38.8 Table 4.1.1 Comparison of modeled values of U /V and Beta, the angle between the surface wind and geostrophic wind, with calcualted values .from Estoque (197.3) and Panofsky and Dutton (1984). 5 0 4 0 --(U 3 0 — B E T A - M O D E L E D v s R N R L Y T I C Zo-2.B0m * o B = l .0 Z o - t . 2 5 m Zo-1.0Gm * 0 Zo°.50m * Zo-.25m & 0 Z o - . 0 8 2 1 m o - ESTOQUE: ( 1 9 7 3 ) * - PANOFSKY and DUTTON ( 1 9 8 4 ) 10 2 0 3 0 40 B E T R R n a l y t i c ( d e g . ) _1 50 I I Figure 4.1.1 F r i c t i o n a l Turning - Modeled vs Calculated -50-large roughness differences over the model domain, the roughness resolution may be as low as five grid spaces or approximately twenty kilometres. This slow response to a step change in roughness is a result of the spatial averaging necessary to calculate the frictional drag at staggered grid nodes. Clearly this limits the ability of the model to respond to the somewhat shorter wavenumber variations found within the real domain. 4.2 Test Run Over a Simple Topographic Feature As stated previously, the one-level model cannot conserve mass. As such, it must rely purely on surface pressure gradients induced by temperature changes to deflect flow around topographic features. To observe this mechanism, the model was run over a very simple topographic feature to see the deflection caused by the feature and to gain insight into the workings of the model. The topographic feature used in the model run was a sinusoidal hill. The cross-section of any radius of the hill is a cosine function with an amplitude of 800m and a wavelength of sixteen grid spacings (see Figure 4.2.1). The value of the roughness length used is .25m and is constant over the domain in order to simplify the situation as much as possible. The upper boundary conditions are such that the slope of the 850mb surface creates a geostrophic flow of 5m/s from 270 degrees near the 850mb level. At the surface upstream from the hill, the flow is slowed and deflected toward the north by frictional effects. As can be seen in Figure 4.2.2, the flow is deflected around the hill. The flow also demonstrates divergence associated with a stagnation point on the windward side of the hill and convergence on the leeward side of the hill. The mechanism behind this behavior of the flow can be seen in Figure 4.2.3. This figure shows the windward dislocation of the surface temperature minimum associated with the peak of the hill. As the flow moves up over the obstruction, potentially cooler air is advected up the side of the hill. This potentially cooler air creates a surface temperature decrease and by so Figure 4.2.1 Sinusoidal H i l l Topography This figure shows the shape of the sinusoidal h i l l used in exploratory modeling. This gri d pattern includes the 3-grid space buffer at the edges of the domain. SINUSOIDAL HILL TIME CHRtMIN) to 0 SURFACE UINDS SCPLE (M/S°)'° w / / '///// • / / / / / ' / //111// / / 1 1 1 / / / . / / t I f / y --- ' / / / / / / / ''/JX/,/,//, A/;M//// i U l h o I Figure 4.2.2 Near Surface Flow Around a Sinusoidal H i l l TEMPERATURE £ TOPOGRAPHY CONTOURS meg. c. e n e t r e s ) INTEGRATION TO STEADY STATE OVER A SINUSOIDAL HILL _i i i i i l_ 4 0 DISTANCE (KH) Figure 4.2.3 Near Surface Temperature Pattern Over a Sinusoidal H i l l . The thin solid lines are temperature contours with a contour interval of 1.0°C. The heavy dashed lines are topographic contours at the 100, 200, 400, 500 and 600 metre elevations. -54-doing, creates a surface pressure minimum. The opposite effect can be observed on the leeward side of the hill where potentially warmer air is advected down the slope. The associated temperature increase causes a local pressure minimum and results in convergent flow behind the hill. This makes the mechanism for the channeling or deflection of flow very much dependent on stability Mass and Dempsey (1985b) point out that with higher velocities and/or lower stabilities the flow will remain relatively undefleeted. Mass and Dempsey suggest that use of the model in high Froude number situations {ie. Fr>1.0) may be inappropriate. This has possible implications when using the model with strong surface heating. Data from Tethersonde flights over the city on several days including the daj'S modeled, show that the potential temperature lapse rate is near 7xlCT3 °K m' 1 . Using the formulation of the Froude number given by Mass and Dempsey (Fr = U/hN), this implies that for the Froude number to be less than 1.0 , U/h must be less than about 10"2s"^. In this formulation, U is the mean wind speed, h is the height of barrier around which the flow is being deflected and N is the Brunt-Vaisalla frequency. If h is chosen to be the height of the mountains surrounding the Fraser Valley, U can be as high as lOm/s and the Froude number will still be less than 1.0. Typically wind speeds are much less than the lOm/s upper limit arrived at here. However, the smaller features on the valley floor are at the same time., of much less vertical extent. Given this upper limit of the Froude number, it seems reasonable to expect the model to produce deflections around the larger valley wall features but to miss the smaller features on the valley floor. This may be of some consequence if validation measurements are made on the valley floor where the model may perform poorly. This problem is examined later in the text. -55-4.3 Modeling Over Idealized Topography The next step in the model testing involved model runs over topography that, in a verj-general way, represented the lower Fraser Valley (see Figure 4.3.1). The intent of this exercise was to determine the effectiveness of the temporal variability of the surface heating parameterization as well as to gain a feel for the model response during integration through an entire day. The surface roughness length (ZQ) was set to lm on the valley bottom, 3.5m over the mountain sides and .01m over the sea surface. The Bowen ratios used were 1.0 for the valley floor, .5 for the mountain sides and .01 for the sea surface. Because of the inexact nature of this aspect of the study, these values were chosen to be near values quoted in the literature and were constant within the three types of surfaces specified. (Cleugh and Oke, 1986; Pielke,1984; Steyn,1980; Oke, 1978; Silversides,1978) The model was integrated to steady state using 850mb pressure and temperature field that were consistent with geostrophic winds near the 850mb level of 5m s~* from 270^. The domain was then cooled using a Q* value of 100W m"w for a period of five hours. There are a number of reasons for initializing the model in this way. The integration of the model to steady state is done adiabatically. Therefore, the steady state fields arrived at are representative of a time during the day when there are no fluxes of heat to or from the atmosphere. The model must be cooled after it reaches steady state in order to bring the fields to a state where they are consistent with those at the beginning of the day being modeled. Therefore the model must be run through part of the previous day. This period of cooling takes place between the time when the surface heat flux goes to zero and the beginning of the day to be modeled. This means that the starting time for the diabatic integration must be a time when heat fluxes at the surface are negligible. Q* reached a value of zero near 19:00 hours on the modeled days. This led to the logical choice of a cooling period of five hours bringing the integration to midnight of the day previous to the one modeled. - 5 6 -Figure 4.3.1 Real and Idealized Topography This figure shows the real (top) and idealized (bottom) discretized topographies used in the modeling study. T h e idealized topography was created by the author to v e r y generally resemble the real topography. -57-The model is cooled using a temporal cosine ramp on the surface cooling. At the beginning of the cooling period, is zero and sinusoidally increases to the 100W m~" value over, thirty time steps. This is done to avoid any step changes in the diabatic forcing and also to represent the smooth change in the surface cooling rate found in reality. Once the domain has been cooled for five hours, the integration continues using the surface heating parameterization described earlier in the text. Presented here are two plots which represent the windfield at 06:00 and 12:00.(see Figures 4.3.2a and 4.3.2b) The 06:00 field shows the land breeze that has developed and also a distinct channeling around the topographic features. There also appears to be a point of stagnation at the mid-point of the foot of the mountain along the northern boundary of the domain. All of these features are what might be expected of a similar situation in reality. The 12:00 frame shows a sea breeze that has set in from the south-west and again the channeling by the large topographic features at the edges of the domain. The rest of the frames through the day, which are not presented here for the sake of brevity, show that the timing of the wind switches over from a land breeze to a sea breeze between the hours of 9:00 and 12:00 . This is consistent with observed times for this switch and also follows the pattern of building from the south-west that is commonly observed (Steyn and Faulkner, 1986). The switch from the sea breeze back to the land breeze occurs somewhat late, in the last three hours of the day. The typical observed times for this switch are nearly three hours earlier (Steyn and Faulkner, 1986). This disparity implies the existence of some thermal inertia over the land that is not modeled or that the specification of the Bowen ratio over the land or sea surface may be in error. This run was made over a very generalized version of the topography and surface characteristics. For this reason speculation of the exact cause of the poor timing of the TIME (HR:MIN) 6: 0 SURFACE WINDS 0.0 10.0 SCALE IM/5) N \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ X \ \ \ I I I / / / I I I / \ \ \ \ \ \ \ I \ \ \ \ \ W \ \ \ — \ \ \ \ / / y / / / / / / / / y / / / / / / / / / / / / / / / / / / / / / / / / / / / / I / / / / / / / / / / / / / / / / / / / / / / / / / / / / / Figure 4 . 3 . 2 a F l o w O v e r I d e a l i z e d T o p o g r a p h y - 0 6 : 0 0 TIME 12: 0 (HR:(1IN) SURFACE WINDS 1 0.0 SCALE 1 10.0 (M/S) y - - - - y ^ • / / / / / / / / / y X / / / / y / / / / / / / y -/ / / / / / / / / y y y y y 1 1 / / / / / / / / / / / / / / / y y y - - — — I I / / I / / / / / / / / / / / / / / / / / y I 1 / / 11 / / / / / / / / / / / / / / / / / / / 1 1 / / 1 i / / / / / / / / / / / / / / / / / / / y I I / / i I / / / / / / / / / / / / / / / / / / / / 1 1 f 1 1 1 / / / / / / / / / / / / / / / / / / / y 1 1 1 j r I / / / / / / / / / / / / / / / / / / / y 1 1 1 1 \ I / / / / / / / / y / / / / / / / / / / y \ \ 1 1 1 1 / / / / / / ^ y y / / 1 1 / / / / / \ \ \ \ \ 1 1 / / / / / y / / I 1 / / / / /' \ \ \ \ 1 1 1 / / / / / — — - / / 1 I / / / / / F i g u r e 4 . 3 . 2 b Flow Over I d e a l i z e d T o p o g r a p h y - 12 :00 - 6 0 -switch, if it occurs with the real topography and surface characteristics, will be held for later when validation methods which are much more diagnostic can be applied. 4.3.1 Edge Filter Effects The use of centred differencing in the numerical scheme used in the model makes it necessary to extrapolate the values of the fields on the perimeter of the domain. The tendencies at the edges of the domain are also set to zero if the flow at the boundarj' is inward. Setting the tendencies to zero is done in order to suppress the propagation of short wavelength perturbations into the modeled fields. It was discovered that this filter applied to the edges of the domain also made the model inappropriately sensitive to the slope of the topography near the edges of the domain. This discovery was made during exploratory modeling over the ideal topography. This filter has the effect of causing spurious deflection of the flow, i.e. up the slopes near the edge of the domain. This occurs where the topography is sloping in such a way as to have a component of the aspect perpendicular to the model boundary. In order to correct this problem, a three grid-space buffer was added to the edges of the domain. The topograph}' of this buffer was horizontal along any line perpendicular to the boundaries of the domain. The addition of this buffer had the desired effect of removing the spurious up-slope winds that were observed without the buffer. Subsequent model runs were all made using the buffer at the edge of the domain. - 6 1 -5. MODELING OVER THE R E A L DOMAIN SURFACE In order to best chose the days to be modeled, some details concerning the synoptic scale features over the domain, as well as their effect on meso-scale flows and circulations must be understood. In the lower Fraser Valley, summer brings an extension of the Pacific anticyclonic regime into the mid-latitudes, resulting in the predominance of high-pressure conditions and generally clear warm weather. Such a pattern may be very persistent, giving prolonged periods of fine weather since any Pacific cyclones are deflected far to the north of the Vancouver area. (Hay and Oke,1976) The synoptic scale pressure gradient over the domain is very weak during these periods of fine weather. This slack synoptic pressure gradient allows the local scale pressure gradients caused by meso-scale variations in surface heating and cooling, to dominate over the synoptic scale background pressure gradient. These local scale pressure gradients result in a number of meso-scale phenomena such as mountain valley winds and land/sea breeze circulations. During the cloudless nights, the radiative balance leads to cooling of the surface and subsequent transfer of sensible heat from the air to the ground surface. This transfer of heat causes a stable layer to form near the surface. This stable stratification within adjoining valleys leads to the creation of mountain valley flows. During clear warm days, land/sea breeze circulations tend to dominate the pattern of flow within the Fraser Valley (Steyn and Faulkner, 1986; McKendry and Steyn, 1986). The flow is topographically constrained to the north and south-east by mountains and channeled to a certain extent in the eastern end of the domain where the valley narrows. Described in a very simplistic fashion, the air is heated over the land surface which in turn creates an upper level pressure gradient from the land to the sea and starts a flow in that direction. Upper level convergence over the sea causes an increase in pressure near the surface and a subsequent sea to land pressure gradient near the surface. At the same time, air over the land surface begins to rise due to instability introduced by surface heating. The near surface pressure gradient causes flow from sea to land to replace the - 6 2 -rising air over the land surface. The result is a circulatory flow with air rising over the land surface and sinking back down over the sea surface. (Atkinson, 1981) This circulation, complete with return flow, has been observed during many soundings done with Tethersonde equipment (Steyn, 1984). Both the mountain valley winds and the land/sea-breeze circulations are profoundly affected by the synoptic scale pressure gradients that exist over the domain (Atkinson, 1981). For example the land/sea breeze effect can appear in more than one form. In the absence of a gradient wind, the land/sea-breeze circulation most often appears in its full circulatory form (Atkinson, 1981). When the synoptic gradient is from the land to the sea, hereafter called an opposing gradient, the sea breeze tends to be more shallow, sets in later and is of shorter duration (Kimble et aZ.,1946; Wexler,1946; Fisher, 1960; Frizzola and Fisher, 1963). In the case where the weak synoptic gradient is from sea to land, hereafter referred to as a following gradient, the sea breeze occurs as an enhancement of the surface level flow from the sea to the land (Sutcliffe,1937). This enhancement of low level wind speed is again the result of the pressure gradient caused by the strong surface heating over the land. This local-scale sea to land pressure gradient is imposed upon the larger scale synoptic gradient at the surface to create this low-level jet within a deep onshore flow. Nocturnal cooling has the effect of retarding the flow at the surface as the local thermally induced surface pressure gradient acts in a direction opposing the synoptic scale gradient. If the synoptic scale pressure gradient is strong enough, it may overpower the local scale pressure gradients completely (Atkinson, 1981). This effect is also exaggerated by the fact that, if the wind speed at the surface is high enough, mechanical mixing will not allow strong temperature gradients to form. In the absence of strong thermal gradients, local-scale pressure gradients cannot form to produce these circulations (Findlater,1963). The pressure gradient caused by larger scale surface heating in the interior of the province has an effect similar to the land/sea heating contrast mentioned above. This - 6 3 -thermal low is formed in very much the same way as the local-scale low surface pressure zone over the lower Fraser Valley, but has effects which are much more wide spread and long lasting. The thermal low of the interior does not have the same intensity of diurnal variation as its analog within the lower Fraser Valley. This thermal low in the interior of the province may dominate the flow pattern over the lower Fraser Valley. As this low is at the surface, its effects at the 850mb upper boundary of the model may not be present. This can cause problems when dealing with a model that only specifies an upper boundary condition which may not account for forcing between the surface and the model's upper boundary. The extent to which the synoptic situation is stationary is also important when using a model of this type. The model assumes that the upper boundary conditions are constant through the duration of the model run. If in reality the 850mb pressure surface or the temperature over that surface changes through the model run, it cannot be expected that the model will simulate the actual surface winds with any accuracy. In order to avoid this problem, modeled days must be chosen that are in the middle of a sequence of consistently fine days so as to minimize the possibility of changes at the upper boundary of the model. Also, in order to assess the model's performance in the two situations, days including both opposing and following synoptic gradients must be chosen. Two such days, July 20,1985 and August 23,1985, were chosen and modeled. What follows is the description of the synoptic situation, the upper boundary conditions used, the model results and the validation of those results for the two days modeled. - 6 4 -5.1 Upper Boundary Conditions The original intent, when choosing the 850mb pressure and temperature surfaces, was to take the Fields from National Weather Centre charts for the region. This approach was tried with very little success. The result of this attempt was to produce winds that were dramatically over-estimated. In retrospect, this should not have been an unexpected result given that the resolution of the charts was nearly two orders of magnitude lower than that of the model as used on the lower Fraser Valley. What proved to give the best results was the estimation of the fields using sounding information taken from pilot balloon flights originating from the Vancouver International Airport. The slope of the 850mb pressure surface is calculated using the measured geostrophic wind and the following formulation for the slope of constant pressure surfaces given by Holton (1979): f V g = k x V eq.5.1 The local lapse rate and the mean temperature over the 850mb surface were calculated using an extrapolation of Tethersonde temperature profiles over the area for that day. One alternative to this is to use information from one of the local upper atmosphere soundings such as Quillayute or Port Hardy. Because these sites are not within the model domain, the data from these sites may vary from the reality of the situation over the domain. It must also be realized that the Tethersonde data will not be available in advance if the model is used in any sort of predictive role, making the use of these alternate upper-air soundings the only choice. The variation of temperature over the 850mb pressure surface was calculated using the local lapse rate and the variation in the 850mb surface elevation over the domain. This gave decidedly small temperature gradients over the upper boundary of the domain, but in the absence of a better estimate this one was used. -65-5.1.1 July 20,1985 As can be seen in the 850mb chart for 1200Z over the Pacific Northwest (Figure 5.1.1), the 850mb pressure surface shows very little slope over the domain. Alhough it has been stated that using these data can lead to drastic over-estimation of 850mb winds, this chart shows the very smooth, low relief of the 850mb surface elevation and temperature fields. The fact that these fields over-estimated geostrophic winds at 850mb would indicate that, if anything, the chart is over-estimating the surface relief. The point to be drawn from this is that the 850mb pressure and temperature fields are very smooth, showing very small gradients over the entire domain. It should also be mentioned that the 700mb chart exhibits these same characteristics. The measured wind at approximately the 850mb level were at 4m/s from 270°. In fact the winds from the surface up to about 2000m only vary 20° to the north or south of due west and, as expected, give very little, if any, indication of baroclinicity. These conditions provided a following gradient necessary for model simulations. 5.1.2 August 23,1985 Again the 1200Z 850mb chart (Figure 5.1.2)shows very little relief in the pressure surface or variation in temperature at 850mb, over the Pacific Northwest. The measured wind at the 850mb level is from 137° at about 4m/s. The wind veers slightly with height above the surface but shearing is nearly nonexistent up to about the 1700m level where the data ends. Although these conditions did not provide a synoptic pressure gradient that directly opposed the local thermally induced daytime pressure gradient, a substantial component of the synoptic gradient was in the east to west direction. These conditions provided the opposing synoptic gradient desired for modeling the thermally driven mechanically forced flows over the domain. - 6 6 -Figure 5.1.1 850mb Chart - 1200Z July 20,1985 This figure shows the elevation in tens of metres (solid lines) of the 850mb pressure surface and the temperature in degrees Celsius (dashed lines) over the 850mb pressure surface. - 6 7 -Figure 5.1.2 850mb Chart - 1200Z August 23,1985 This figure shows the elevation in tens of metres (solid lines) of the 850mb pressure surface and the temperature in degrees Celsius (dashed lines) over the 850mb pressure surface. -68-5.2 Statistics Used in Model Validation There are a number of statistical indices that can be used to assess the performance of a model. Willmott, (1981) in a paper on the validation of models, recommends using a select number of these indices for model validation. Willmott (1981) recommends the calculation of observed and predicted means and standard deviations as well as the slope and intercept of a least-squares regression between the predicted and observed values. Also included in the recommended list are the root mean square error (RMSE), the systematic root mean square error (RMSES), the unsystematic root mean square error (RMSEU) and the index of agreement. With the exception of the slope and intercept of the least squares regression between predicted and observed values, this list of indices has been adopted. The formulations for these terms is presented in Appendix VI. With one exception, this list of indices is comprised of standard statistical terms. The exception is the index of agreement. Willmott (1984) finds r, the correlation coefficient and 2 r to be inappropriate statistics, at best, for making quantitative comparisons between observed and predicted fields. Willmott et al. (1985) make a very sound case for the use of the index of agreement in the place of the coefficient of determination for model validation. Making use of this reasoning, the index of agreement is used in place of the coefficient of determination as a summary index of model performance. The above statistics are calculated using modeled values which are inverse square, weighted averages of the U and V components of the four grid nodes surrounding any observation site. The U and V components are treated statistically as separate independent estimates. This is can be done because the model calculates the U and V components separately, although the U and V components are not totally independent of each other. - 6 9 -5.3 Modeled and Observed Fields and Validation Data On both the chosen days the model was integrated to steady state using the upper boundary conditions specified. The domain was subsequently diabatically cooled for five hours in order to bring the fields to a state consistent with the beginning of the modeling day. The integration was then allowed to proceed for a period of 24 hours using the diabatic heating and cooling specified by the modeled thermal regime of that particular day. The results are presented here in the form of hourly averaged windfields (modeled and observed), as well as the wind residual plot and statistical indices for that hour. The wind vector residual plot is a diagnostic tool that allows the visual assessment of model performance. The observed wind vector is subtracted from the modeled wind vector at all of the validation points. On the wind vector residual plots, the station number with a circle around it is plotted at the head of each wind velocity residual vector. The tails of the vectors are at the centre of the map and the concentric circles indicate the magnitude of the residual wind vector. The purpose of this is to point out any systematic errors in the modeled field and to show the degree to which the residual vectors are scattered about the perfect prediction at the centre of the diagram. (See Figures 5.4.1-24c and 5.5.1-24c) -70-5.4 Discussion of Modeling of July 20,1985 As can be seen in Figure 5.4.1a, after the first hour of integration, the model produces winds that are southerly at the mid-point of the southern boundary of the domain. The flow diverges to the west and east at points north of the southern boundary. The flow tends to stagnate against the northern mountains. This stagnation is to be expected given the stable stratification of the atmosphere in the lowest 1000m. The flow does show some channeling in the north-west corner of the domain where it is forced northward into Howe Sound. The offshore flow in the western half of the domain is being caused by the differential heating of the land and sea surfaces. In the eastern half of the domain, the flow appears to be dominated by the synoptic gradient. The observed winds show that the actual flow at this time is from the east in the eastern-most part of the domain, (see Figure 5.4.1b) It appears that this discrepancy is, at least in part, a result of the fact that the local pressure differential created by the differential heating has effects that do not reach to the eastern border of the domain. There are two main reasons why the effects only reach to just past the mid-way point of the domain. The explanation of these possibilities will serve to elucidate two characteristics of numerical models of this type. The role of temperature diffusion in a numerical model of this nature is to simulate sub-grid scale fluxes such as turbulent transport of heat and to suppress numerical instability (Mass and Dempsey, 1985b; Pielke,1984). The amount of diffusion specified to do the job of suppression of numerical instability may be somewhat larger than its real counterpart, which only represents the sub-grid scale fluxes. Also, the coefficient of diffusion is specified as constant over the entire domain when, in reality, it may be inhomogeneous and anisotropic (Panofsky and Dutton, 1984). The end result is the specification of a coefficient of diffusion which may be too large. This will have the effect of "smearing" the heat over the domain, thus reducing temperature contrasts. - 7 1 -The second of the two mechanisms that may be responsible for the discrepancy can be found in the fact that the model has no way of estimating upstream conditions. Inflow at the boundaries tends to flood the domain with air that has a potential temperature equal to that at the boundary. In this case, the flow comes in from the southern boundary and tends to eradicate the temperature differential built up by differential diabatic cooling, (see Figure 5.4.1-4a) As the model is integrated from steady state, the differential land/sea cooling creates a pressure gradient across the coastline causing the flow to slowly swing around from a westerly to a southerly. The balance of pressure gradient, friction and Coriolis forces, dictates a strong southerly flow that effectively isolates the inland portion of the domain from the land/sea temperature contrast, (see Figure 5.4.1a) This flood of air at a constant potential temperature is only affected by the temperature differential close to the shoreline where the contrast is still being created by the differential cooling. Further inland, the flow is driven by the synoptic scale pressure gradient in the absence of a local, thermally produced gradient. In the real domain, the coastline swings eastward south of the domain. This implies that there would in reality, be a component to the local thermally created pressure gradient that might oppose the rapid inflow at the southern boundary. This possibility seems to be borne out by the fact that the flow at the mid-point of the southern boundary is observed to be very light, (see Figure 5.4.1b) This points to the modeled strong southerly flow being an artifact of edge effects. Observed flow at the eastern boundary is more likely katabatic outflow at the surface, possibly separated from an upper level flow occurring in the opposite direction, (see Figure 5.4.1b) Again, because the model has no way of estimating upstream conditions, this outflow is not represented by the model. These two factors serve to illustrate some of the problems encountered when conditions outside the boundary of the model affect the flow inside the domain. One - 7 2 -possible solution to this problem might be to telescope the model grid away from the present .boundaries. This could be done by adding another few rows around the domain with a much larger grid spacing in order to better estimate upstream conditions. The plot of the wind velocity residuals shows a consistent over-prediction of the southerly component of the flow and a low index of agreement, (see Figure 5.4.1c) This rather poor model performance is a result of the above mentioned limitations of the model used over this particularly complex domain. The modeled and observed flow patterns remain very much unchanged for the next four hours, (see Figures 5.4.2-5a) The model performance remains poor but improves slowly. At the end of the sixth hour of simulation, the outflow from the surrounding valleys has stopped and the synoptic gradient has become dominant in the far eastern portion of the domain, (see Figure 5.4.6a) The model is still over-predicting the southerly component of the flow but the agreement in the eastern half of the domain is dramatically improved. The index of agreement increases accordingly to a value of .5598. By the end of the seventh hour of simulation (see Figure 5.4.7a) the radiative cooling at the surface has given way to positive values of Q* and the temperature gradient across the coastline is beginning to dissipate. This results in the reduction of the land-to-sea pressure gradient and a clockwise rotation of the flow in the western half of the domain. This rotation of the flow is predicted by analytical studies of the sea breeze presented by Haurwitz (1947). Further analytical, numerical and observational studies by Neumann (1977), Hsu (1970) and Angell and Pack (1965) show this same rotation. This modeled rotation of the flow continues as the surface heating increases through the morning. As the flow rotates, the agreement increases to its maximum value of .8339 after the ninth hour of simulation. By the eleventh hour of integration (see Figure 5.4.11a), the model produces a predominantly south-westerly flow in the western half of the domain as the sea breeze slowly fills in. - 7 3 -The effects of the complex coastline can be seen in the flow pattern over Boundary Bay, immediately south of Vancouver, (see Figure 5.4.11a) Here the flow appears to be divergent over the bay. McPherson (1970), in a numerical study of the effects of a coastal irregularity on the sea breeze, shows modeled flows that exhibit this same divergence over a hypothetical bay. Although the vertical motion implied by this divergence at the surface cannot be verified, this modeled feature serves as an indication that the model can resolve spatial variations in surface heating on the scale of at least the size of Boundary Bay. Temporally, the model produces the onset of the sea breeze earlier than the observed onset. By noon the model predicts a nearly fully developed sea breeze when in reality, this does not fully set in until around 13:00. This premature, clockwise swing of the flow draws the agreement down to a value of approximately .45, where it stays until nearly 15:00. Once the modeled sea breeze has set in, the flow appears to be close to steady state with the speed of the flow decreasing near the eastern boundary. Observational studies over flat topography in the mid-latitudes show that as the sea breeze grows in intensity throughout the day, it swings to a direction parallel to the coast (Atkinson, 1981). The large westerly component in the observed flow appears to be, at least in part, the result of deflection created by the mountains on the northern boundary of the domain. The model seems to do a good job of reproducing this channeling effect. The observed flow near the eastern boundary has a much larger speed than the modeled flow. The reason for the observed flow being as large as it is, may be found in a larger scale phenomenon. As previously mentioned, the larger scale thermal low in the interior of the province may be drawing the near-surface flow up the valleys toward the interior of the province. If this is the case, then it demonstrates another good reason for telescoping the boundaries to include a rough representation of such large scale effects. Figure 5.4.25 shows the index of agreement through the modeled July day. The indices of agreement are for the run being discussed (spatially varying surface characteristics), a run with spatially constant surface characteristics and a run with attenuated gradients in - 7 4 -the upper boundary conditions. As can be seen upon examination of curve number 1, the agreement starts to drop at about 15:00 to reach a minimum at 19:00 and then begins to rise to values near .50 at the close of the day. This dip in the level of agreement is caused by the inability of the model to reproduce the timing of the decay of the sea breeze and the onset of the land breeze. The westerly component of the modeled flow only seems to start to decay at about 19:00. (see Figure 5.4.19a) The observed flow shows signs of reversal in the extreme western portion of the domain as early as 18:00. (see Figure 5.4.18b) The westerly component of the observed flow decays until near 21:00, when the flow has an easterly component over most of the domain, (see Figure 5.4.21b) Although the westerly component of the flow over the eastern half of the domain has decreased substantially, the modeled flow is still predominantly westerly. It is not until the last hour of simulation that the flow begins a more rapid anticlockwise swing around to the west, (see Figure 5.4.24a) This lag in the modeled transition from a sea to land breeze can be seen in the vector residual plots, (see Figures 5.4.17-22c) At 17:00 the vector residuals are scattered about the centre of the plot in a somewhat uniformly distributed fashion. As the integration progresses through the afternoon, the residual plots tend to be distributed to the right of centre. At 22:00 all of the vectors are on the right side of the plot. This is a clear indication that the model is over-predicting the westerly component of the flow. The early onset of the sea breeze and the late reversal that the model produces may be caused by improper specification of the upper boundary conditions. Investigation of this possibility can be carried out by examining model response to changed upper boundary conditions. If the synoptic gradient were attenuated, the local scale pressure gradient may not have to be as strong to reverse the flow in the evening and the enhancement of the westerly flow in the mid-morning might be diminished. Figure 5.4.25 shows the index of agreement for a model run using an 850mb pressure surface with the same shape as the original run but with the 850mb surface slope and the -75-850mb temperature gradient attenuated by a factor of 10. (see Figure 5.4.25,curve 2) As is shown, the results were not significantly different. This indicates that the local pressure gradients produced by differential heating and cooling are the dominant forcing in the modeled windfields. This dominance of local heating and cooling points to the possibility that the modeled heating and cooling in the eastern half of the domain may be under-estimated. The results of a model run with the same upper boundary conditions as the first run, but with constant surface characteristics, prove interesting. In this run, a roughness of ZQ=1.0m and a Bowen ratio of 1.0 are specified over the entire land surface. As can be seen by comparing this to Figure 3.1.1, this leaves the Bowen ratio nearly unchanged over the valley floor. The significant change in the Bowen ratio is over the mountain slopes where it is increased to a value of 1.0. This is an increase of the Bowen ratio by a factor of 2.0 over a substantial part of the domain. Curve number 3 of Figure 5.4.25 shows the index of agreement through the day for this run. The agreement does not show the same strong peak near 08:00 exhibited by the run using spatially varying surface characteristics. There is however, a significant increase in agreement near the time when the observed flow is switching over from a sea breeze to a land breeze. The plotted windfields for this period show that the modeled flow reverses much earlier, which is more in agreement with the observed flow pattern. The reason for this would appear to be the increased diabatic forcing over the mountain slopes as compared to the initial runs. It appears that the enhanced cooling over the mountain slopes is enough to reverse the flow at the proper time, improving the agreement at the time of switch-over from the sea to land breeze. This change is not without some disadvantages in that it also has the effect of decreasing the agreement in the morning hours, (see Figure 5.4.25) These experiments indicate a number of things. The first and probably most important is that the model is very sensitive to spatial variability of surface characteristics, at least - 7 6 -in some synoptic situations. This coupled with the fact that the agreement does not increase when the synoptic gradient is attenuated indicates that the lack of agreement may well be caused by the model's response to surface heating rather than the improper-specification of the synoptic gradient. This is encouraging from the point of view of fine-tuning some of the systematic error out of the model. A number of experimental model runs under varied synoptic conditions could be made in order to gain a better understanding of the model response to specific spatial patterns of surface characteristics. Hopefully, the improved specification of the spatially varying Bowen ratio would lead to improved agreement. It is also worth noting that if the model is ever to be used in a predictive mode, the model operator will not have the advantage of real time data to compare the modeled fields with. In all likelihood, the only available data may be the upper-level geostrophic winds and a temperature profile from outside the domain. The severity of this problem may be minimized by the fact that in situations where a slack synoptic gradient exists over the model domain, the model may be fairly insensitive to errors in upper boundary condition specification. -77-Figures 5.4.1-24a These figures show the modeled windfields over the domain for the 24 modeled hours of July 20,1985. The end of the one hour averaging periods are shown in the top left corner of the frames. The wind speed scales appear in the top right corner of the plot frames. Figures 5.4.1-24b These figures show the interpolated windfields for the 24 modeled hours of July 20,1985. The windfields shown are inverse distance weighted least squares fits of 2nd order polynomials to the data points. Only the sections of the windfields that are close to data collection sites are shown. This is because the reliability of the interpolation decreases rapidly with distance from measured points. Figures 5.4.1-24c These figures show the windfield residual plots for the 24 modeled hours of July 20,1985. The wind velocity residuals are the result of subtracting the observed from the predicted wind vectors. Circles enclosing the site numbers are plotted at the tip of the residual vectors with the tails of the residual vectors at the centres of the concentric circles. -78-Figure 5.4.1a Figure 5.4.1b - 7 9 -TIME. HOUR ENOING 1:00 RMSE d / s ) = 1.8192 RMSES ( i / s ) = 0.8837 RflSEU (n/s) = 1.5901 000 deg. 270 deg. D2 = 0.3832 MEAN OBS. ( i / s ) = 0.491 MEAN PREO. d / s ) = 1.51 SIGMA OBS. = 0.7498 SIGMA PRED. = 1.6271 090 deg ~ i 1 0.0 3.0 6.0 SCALE (M/S) 180 deg. WIND VbLUCIlY RhSlUUALS Figure 5.4.1c TIME. HOUR ENDING 2:00 RMSE (m/s) = 1.7077 RMSES («/s) = 0.7062 RMSEU («/s) = 1.5548 000 deg. 270 deg. 02 = 0.2868 MEAN OBS. («/s) = 0.337 MEAN PREO. ( i / s ) = 1.49 SIGMA OBS. = 0.4910 SIGMA PRED. = 1.5926 090 deg. 0.0 3.0 6.0 SCALE (M/S) 180 deg. WIND VbLUCIlY RbSlUUflbS Figure 5.4.2c Figure 5.4.2c - 8 0 -JUL_20_B5 Ur£ (HRtMIN) 2i 0 S U R F A C E W I N D S o.o ip .o SCPLE (H/S) « so eo OISTRNCE (KM) Figure 5.4.2a 1985.07.20: 200 i i i Figure 5.4.2b - 8 1 -J U L _ 2 0 _ 8 5 T T r E (HRtfUN) 3i 0 S U R F A C E W I N D S 10O0M CONTOUR (M/ 10 .a S) a s-—T" 10 » \ \ \ \ \ \ \ \ \ \ \ \ / s / s / y / s / s / ' / / / / <7 20 4 SO 80 OISTRHCE (KM) ao n i n F i g u r e 5 . 4 . 3 a 1985.07.20: 300 i i , i F i g u r e 5 . 4 . 3 b - 8 2 -TIHE. HOUR ENDING 3:00 RMSE d/s) = 1.8013 RMSES d/s) = 0.9083 RMSEU (n/s) = 1.5556 000 deg. 270 deg. 02 = 0.3985 MEAN OBS. (n/s) =• 0.490 MEAN PRED. (i /s) = 1.51 SIGMA OBS. = 0.8187 SIGMA PRED. = 1.6098 090 deg. - i 1 0 . 0 3 .0 6 .0 SCALE (M/S) 180 deg. UIND VELOCITY RESIDURLS Figure 5.4.3c TIME. HOUR ENDING 4:00 RMSE (n/s) = 1.7936 RMSES (m/s) = 0.8652 RMSEU (m/s) = 1.5711 000 deg. 270 deg. 02 = 0.3539 MEAN OBS. (a/s) = 0.405 MEAN PRED. (i /s) = 1.53 SIGMA OBS. = 0.7024 SIGMA PRED. = 1.6275 090 deg. T 0 .0 3 .0 SCALE (M/S) 180 deg. i 6.0 UIND VELOCITY RESIDUALS Figure 5.4.4c Figure 5.4.4c - 8 3 -JUL_20_85 TIME (HRtHIN) 4i 0 S U R F A C E W I N D S lOOOtl CONTOUR 10.0 S) 40 SO 80 OISTfiNCE (KM) Figure 5.4.4a 1985.07.20: 400 Figure 5.4.4b -84-JULY_20J98S TIME (HRifUN) 5 i 0 S U R F A C E U I N D S 1000H CONTOUR o.o 10.0 SCflLE IK/SI s » f / / 1 1 / / / / / 1 / / / / S s 1 / / / / / 1 / / / / s * t / / / / y 1 / / / / s / f / / / / / / [ / / / / / l 1 60 70 80 I / / I \ \ *1 \ DISTANCE (KM) too Figure 5.4.5a 1985.07.20: 500 Figure 5.4.5b -85-TIME. HOUR ENDING 5:00 RMSE d/s) = 1.5367 RMSES (i/s) = 0.6157 RMSEU («/s) = 1.4079 000 deg. 270 deg. 02 = 0.4793 MEAN OBS. (»/s) = 0.392 HERN PRED. («/s) = 1.52 SIGMR OBS. = 0.6103 SIGMR PRED. = 1.6296 090 deg. T 1 0.0 3.0 6.0 SCALE (M/S) 180 deg. UIND VELOCITY RESIDUALS Figure 5.4.5c TIME. HOUR ENDING 6:00 RMSE («/s) = 1.5267 RMSES (d/s) = 0.6734 RMSEU (»/s) = 1.3702 000 deg. 270 deg. D2 = 0.5598 MEAN OBS. d/s) = 0.512 MERN PRED. (i/s) = 1.54 SIGMR OBS. = 0.7575 SIGMR PRED. = 1.6421 090 deg. 0.0 3.0 6.0 SCALE (M/S) 180 deg. iJIND VELOCITY RESIDUALS Figure 5.4.6c - 8 6 -JULY_20.1985 TinE (HRiMIN) 61 0 SURFACE WINDS 10ODM CONTOUR 1 1 0.0 10 0 SCRLE (M/S) 1 r *" OISTRNCE (KM) F i g u r e 5 . 4 . 6 a 1985.07.20: 600 .,, , 1... ' . i . F i g u r e 5 . 4 . 6 b -87-J U L Y _ 2 0 _ I 9 8 S TIME (HRiMIN) It 0 , S U R F A C E U I N D S 0.0 10.0 SCALE (M/S) • /TV* Figure 5.4.7a Figure 5.4.7b - 8 8 -TIME. HOUR ENDING 7:00 RMSE d/s) - 1.3826 RMSES («/s) = 0.8387 RMSEU («/s) = 1.0992 000 deg. 270 deg. 02 = 0.6774 MEAN OBS. l»/s) = 0.609 MEAN PRED. («/s) =• 1.52 SIGMA OBS. = 0.8483 SIGMA PRED. = 1.6104 090 deg. ~i 1 0.0 3.0 6.0 SCALE (M/S) 180 deg. UIND VELOCITY RESIDUALS F i g u r e 5 . 4 . 7 c TIME. HOUR ENDING 8:00 RMSE (i/s) = 1.3054 RMSES (•/s) = 0.8051 RMSEU im/s) = 1.0276 000 deg. 270 deg. D2 = 0.7707 MEAN OBS. ( i /s ) = 0.866 MEAN PRED. (n/s) = 1.59 SIGMA OBS. = 1.0952 SIGMA PRED. = 1.6648 090 deg. n 1 0.0 3.0 6.0 SCALE (M/S) 180 deg. UIND VELOCITY RESIDUALS F i g u r e 5 . 4 . 8 c -89-JULY_20_I985 TIME (HRiMIN) 81 0 SURFACE UINDS 1000M CONTOUR o.o ip.o SCflLE (M/S) 40 SO 60 DISTANCE (KM) Figure 5.4.8a Figure 5.4.8b -90-J U L Y _ 2 0 _ 1 9 8 5 TIME (HRiHIN) 9 i 0 SURFACE UINDS 1000M CONTOUR SCflLE tM/s") 40 SO 80 OISTfiHCE (KM) Figure 5.4.9a 1985.07.20: 900 Figure 5.4.9b - 9 1 -TIME. HOUR ENDING 9:00 RUSE d/s) - 1.0429 RMSES («/s) =• 0.5900 RMSEU l»/s) = 0.8599 000 deg. 270 deg. D2 - 0.8339 MERN OBS. («/s) - 1.192 MEAN PRED. (i /s) = 1.36 SIGMR OBS. = 1.3372 SIGMP PRED. = 1.3170 030 deg. i 0.0 3.0 6.0 SCALE (M/S) 180 deg. JIND VELOCITY RESIDUALS Figure 5.4.9c TIME. HOUR ENDING 10:00 RMSE ta/s) - 1.5925 RMSES d/s) = 1.3733 RMSEU (>/s) - 0.8062 000 deg. 270 deg. D2 = 0.5801 MEAN OBS. d/s) = 1.233 MERN PRED. (PVS) = 1.35 SIGMR OBS. = 1.5426 SIGMR PRED. = 0.9221 090 deg. 0.0 3.0 6.0 SCRLE (M/S) 180 deg. UiIND VELOCITY RESIDUALS Figure 5.4.10c - 9 2 -JULY_20_1985 TIME (HRIMIN) 10: 0 S U R F A C E U I N D S 10ODH CONTOUR SCALE (M/s")' o \ \ \ \ \ \ \ \ \ i \ \ \ \ \ \ w w w w\\\\ T 40 30 10 DISTANCE (KM) Figure 5.4.10a 1985.07.20:1000 i i i Figure 5.4.10b - 9 3 -JULY_20_1985 TIME (HRtMINl lit 0 S U R F A C E W I N D S 1000M CONTOUR SCflLE (M/S) it.a -r « S3 90 OISTRHCE (KM) F i g u r e 5 . 4 . 1 1 a 1985.07.20:11.00 F i g u r e 5 . 4 . 1 1 b -94-TIME. HOUR ENDING 11:00 RUSE («/s) - 1.9360 RMSES d/s) =• 1.7699 RMSEU («/s) = 0.7845 000 deg. 270 deg. D2 = 0.4874 MEAN OBS. ( i/s) = 1.322 MEAN PRED. ( i/s) = 1.75 SIGMA OBS. = 1.5054 SIGMA PRED. =» Q.8150 090 deg. 0.0 3.0 6.0 SCALE IM/S) 180 deg. LJIND VELOCITY RESIDUALS Figure 5.4.11c TIME. HOUR ENDING 12:00 RMSE («/s) = 1.8328 RMSES («/s) = 1.6278 RMSEU (a/s) = 0.8421 270 deg. 02 = 0.4414 MEAN OBS. (n/s) - 1.440 MEAN PRED. ( i/s) = 2.01 SIGMA 08S. •- 1.3908 SIGMA PRED. = 0.8617 000 deg. 090 deg. 0.0 3.0 6.0 SCALE (M/S) 180 deg. UIND VELOCITY RESIDUALS Figure 5.4.12c -95-JULY_20_1985 TIME (HR.HIN) 12i 0 , S U R F A C E U I N D S 1000M CONTOUR 10.0 £ (M/S) ' / / / / / / y / / s a-/ / ^ 6 7 / / // / / / / / ' ' ' 7'/ / J r i / ' ' ' ' ' 7 20 49 SO 80 DISTANCE (KM) ao Figure 5.4.12a 1985.07.20:1200 IOOOH CONTOUR &flLE (M/S*) * SBTFLED PtttKTS * OISTM& (toil Figure 5.4.12b -96-J U L Y . 2 0 _ 1 9 8 S TIME (HRihTN) 13< 0 S U R F A C E W I N D S lOOOtl CONTOUR (M/S) ° / / %''l7... 40 S3 SD 0ISTONCE (KM) Figure 5.4.13a Figure 5.4.13b -97-TIME. HOUR ENOING 13:00 RMSE («/s) = 1.5361 RMSES (m/s) = 1.2599 RMSEU («/s) = 0.8787 000 deg. 270 d e g . D2 = 0.4562 MEAN OBS. («/s) - 1.608 MEAN PRED. ui/s) = 2.03 SIGMA OBS. = 1.2161 SIGMA PRED. = 0.9031 090 deg. T 0.0 3.0 6.0 SCALE (M/S) 180 deg. UIND VELOCITY RESIDUALS Figure 5.4.13c TIME. HOUR ENDING 14:00 RMSE («/s) = 1.4953 RMSES («/s) = 1.1949 RMSEU (»/s) = 0.8990 000 deg. 270 d e g . D2 = 0.4637 MEAN OBS. (i /s) - 1.753 MEAN PREO. (i /s) • 1.96 SIGMA OBS. = 1.2917 SIGMA PRED. = 0.9265 090 deg. ~r 1 0.0 3.0 6.0 SCALE (M/S) 180 deg. UIND VELOCITY RESIDUALS Figure 5.4.14c -98-JULY_20_1985 Tire, mitnm S U R F A C E UINDS 1000M CONTOUR SCflLE 10.0 (M /S) / / / / L /// ///'/ / / / / / DISTANCE (KM) " Figure 5.4.14a 1985.07.20:1400 * OISTOwk (KM) " Figure 5.4.14b - 9 9 -.20.1985 (HR.MINJ SURFACE UINDS 10O0M CONTOUR 10.0 (K/S) / / / / / / / / / / / z / / / / / / / / / / / / WW///. 1 4) so n DISTANCE (KM) Figure 5.4.15a 1985.07.20:1500 Figure 5.4.15b -100-TTME. HOUR ENDING 15:00 RMSE (m/s) = 1.4710 RMSES im/s) = 1.1885 RMSEU (n/s) = 0.8668 270 deg. D2 = 0.5249 MEAN OBS. (i /s) = 2.065 MEAN PRED. im/s) = 1.91 SIGMA OBS. = 1.4366 SIGMA PRED. = 0.9107 000 deg. 090 deg. i 0.0 3.0 6.0 SCALE (M/S) 180 deg. WIND VELOCITY RESIDUALS Figure 5.4.15c TIME. HOUR ENDING 16:00 RMSE (n/s) = 1.6449 RMSES (»/s) = 1.4065 RMSEU (a/s) = 0.8531 270 deg. D2 = 0.3802 MEAN OBS. im/s) = 2.067 MEAN PRED. (n/s) = 1.87 SIGMA OBS. = 1.4983 SIGMA PRED. = 0.8695 000 deg. 090 deg. -i 1 0.0 3.0 6.0 SCALE (M/S) 180 deg. WIND VELOCITY RESIDUALS Figure 5.4.16c -101-JULY_2Q_1985 TIME (HRiMIN) ISi 0 SURFACE UINDS 0.0 10.0 SCALE (M/5) 'IS?)} ^ / / ss// ss///A ////// //////. Www-DISTANCE (KM) " Figure 5.4.16a 1985.07.20:1600 , i , , i i Figure 5.4.16b -102-J U L Y _ 2 0 _ I 9 8 5 TIME (HRiMIN) 17: 0 SURFACE UINDS I000M CONTOUR SCALE (M/S0)'0 S a / / / / /// //// 4 so n DISTANCE (KM) Figure 5.4.17a 1985.07.20:1700 Figure 5.4.17b -103-TIME. HOUR ENDING 17:00 RMSE («/s) = 1.6753 RMSES («/s) = 1.4695 RMSEU Im/s) = 0.8045 000 deg. 270 deg. 0 2 = 0.3092 MEAN OBS. Cm/s) » 1.823 MEAN PRED. (i/s) = 1.85 SIGMA OBS. = 1.3715 SIGMA PRED. = 0.8235 090 deg. 0.0 3.0 6.0 SCALE (M/S1 180 deg. WIND VELOCITY RESIDUALS Figure 5.4.17c TIME. HOUR ENDING 18:00 RMSE Im/s) = 1.6797 RMSES («/s) = 1.5027 RMSEU d/s) = 0.7506 000 deg. 270 deg. D2 = 0 .2924 MEAN OBS. (n/s) = 1.501 MEAN PRED. (i/s) - 1.79 SIGMA OBS. = 1.2737 SIGMA PRED. = 0.7771 090 deg. T 1 0.0 3.0 8.0 SCALE (M/S) 180 deg. u/IND VELOCITY RESIDUALS Figure 5.4.18c -104-JULY_20_1985 TIME (HRifUN) 18! 0 SURFACE WINDS 10O0H CONTOUR o.g lo.o SCALE (M/S) / / / / ss / / / ////// ////// i 40 so so DISTANCE (KM) Figure 5.4.18a 1985.07.20:1800 i i i Figure 5.4.18b -105-SURFRCE WINDS / / / / / / y / / / / y y y y / / / y y ~~ 20 « 50 SO DISTPNCE (KM) ao 90 Figure 5.4.19a 1985.07.20:1900 11X01 CONTOUR icns (n/s*) * sfm£D POINTS * DISTANCE HOI) Figure 5.4.19b - 1 0 6 -TIME. HOUR ENDING 19:00 RMSE d/s) = 1.6238 RMSES d/s) = 1.5055 RMSEU (»7s) - 0.6084 000 deg. 270 deg. 02 = 0.2931 MERN OBS. (i/s) = 1.244 MEAN PRED. (i/s) = 1.51 SIGMR OBS. = 1.2104 SIGMA PRED. = 0.6448 090 deg. ~\ 1 0.0 3.0 8.0 SCALE (M/S) 180 deg. WIND VELOCITY RESIDUALS F i g u r e 5 . 4 . 1 9 c TIME. HOUR ENDING 20:00 RMSE (>/s) = 1.5780 RMSES (i/s) = 1.5163 RMSEU (n/s) = 0.4369 000 deg. 270 deg. 02 = 0.3536 MEAN OBS. (»/s) = 1.012 MERN PRED. (i/s) = 1.21 SIGMA OBS. = 1.1680 SIGMR PRED. = 0.4644 090 deg. 0.0 3.0 6.0 SCALE (M/S) 180 deg. WIND VELOCITY RESIDUALS F i g u r e 5 . 4 . 2 0 c -107-J U L Y _ 2 0 _ I 9 8 5 TIME (HRiMIN) 20. 0 , SURFACE UINDS i 1 0.0 10.0 SCALE (M/S) 2*i o ^ * s / / / "V / ^ ^ y y y / / 7/ y / / / / / / / / / A / / / / / / / V / / / / / / / / / / / / / / / / / / / / / U - / 40 50 80 DISTANCE (KM) Figure 5.4.20a 1985.07.20:2000 Figure 5.4.20b -108-JULY. TIME 21i 0 .20J98S (HRtMIN) SURFACE UINDS y y y / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / „ t t t t l t / / / * C " ' " - i r i 40 so n DISTANCE (KM) Figure 5.4.21a 1985.07.20:2100 Figure 5.4.21b - 1 0 9 -TIME, HOUR ENOING 21:00 RUSE («/s) = 1.6526 RMSES (i/s) » 1.6137 RMSEU (m/s) = 0.3563 000 deg. 270 deg. 02 = 0.3643 MEAN OBS. (m/s) = 0.620 MEAN PREO. (m/s) = 1.07 SIGMA OBS. = 0.9622 SIGMA PRED. = 0.3765 090 deg. 0.0 3.0 6.0 SCALE (M/S) 180 deg. UIND VELOCITY RESIDURLS Figure 5.4.21c TIME. HOUR ENDING 22:00 RMSE (m/s) = 1.3591 RMSES (m/s) = 1.3017 RMSEU (m/s) = 0.3908 000 deg. 270 deg. D2 = 0.4139 MEAN OBS. (i/s) = 0.474 MEAN PRED. (i/s) = 1.03 SIGMA OBS. = 0.7042 SIGMA PRED. - 0.3984 090 deg. 0.0 3.0 6.0 SCALE (M/S) 180 deg. UIND VELOCITY RESIDUALS Figure 5.4.22c MX) 30NUI5I0 OB OS 0> / / 1 / / / I / / / J / / / J / / / / / / / / / / / / / / / / / / / y / / y y / </ / y / y / / 7, fi > o 2 o m -H 4- 1 ^ ^ (S/H) 3THDS ooi T o anoiNCO uoooi SONIfl 33bdcJnS (NIW'dH) -on-— I l l — JULY_20_1385 TIME (HRittlN) 23i 0 SURFACE UINDS 100O1 CONTOUR 0.0 10.0 SCfllE (M/S) T 40 SO 80 DISTANCE (KM) Figure 5.4.23a 1985.07.20:2300 III I ' U . ' 1 I ' .11 I i i f i r * DI5TM& (KM! " Figure 5.4.23b -112-TIME. HOUR ENDING 23:00 RMSE d/s) - 0.9998 RMSES («/s) = 0.9015 RMSEU d/s) = 0.4322 000 drg. 270 deg. D2 = 0.5115 MEAN OBS. (n/s) = 0.413 MEAN PRED. (n/s) = 1.03 SIGMA OBS. = 0.6163 SIGMA PRED. = 0.4789 090 deg. i 1 0.0 3.0 8.0 SCALE (M/S) 180 deg. WIND VELOCITY RESIDUALS Figure 5.4.23c TIME. HOUR ENDING 24:00 RMSE (n/s) = 1.4000 RMSES (II/S) = 1.2761 RMSEU (n/s) = 0.5758 000 deg. 270 deg. D2 - 0.4465 MEAN OBS. (i/s) - 0.555 MEAN PRED. (n/s) = 1.01 SIGMA OBS. = 1.0011 SIGMA PRED. = 0.5918 090 deg. i 1 0.0 3.0 6.0 SCALE (M/S) 180 deg. WIND VELOCITY RESIDUALS Figure 5.4.24c -113-SURFRCE UINDS wan CONTOUR o.g ip.o SCALE (M/S) / / / / / / / / / / / / / / / / UL / / / / / / / / / / / / / / / / / / / / / / / M 10 40 SO 60 DISTANCE (KM) Figure 5.4.24a 1985.07.20:2400 SCTU IM/?; Figure 5.4.24b INDEX OF RGREEMENT cn cn n C TJ Z X) ID a TJ TJ n I—I "n n co -< T J 2 a T J X) a 3: m I i -r 1 i I I t 4-i T T T - 7 T T -- 1 1 5 -5.5 Discussion of Modeling of August 23.1985 O n August 23, 1985 the synoptic gradient is opposing the land-to-sea pressure gradient created by differential heating. The flow at the surface is f rom the east over most of the domain. In the plots of the windfield in the early morning hours, the effects of the land/sea contrast in surface roughness are apparent, (see Figure 5.5.1a) The flow over the sea surface has a much greater speed and is not turned by the friction as much as the flow over the land surface. It is difficult to see any var iat ion in direction or speed over the land that might be taken as evidence that the varied roughness over the land has any substantial effect. The only exception to this might be the flow paral le l to the mountains near the base of the mountain slopes in the south-east corner of the domain. A l o n g these slopes the modeled roughness is high (ZQ = 3.5m).This high roughness appears to decrease the flow speed along the base of the mountains . This effect m a y in part be caused by variat ions in terrain height but one would expect some noticeable var ia t ion in direction along the base of the slope i f this were the case. A l o n g the base of the slopes the flow is very smooth. Channel ing can again be observed i n the eastern half of the domain where the flow diverges as it flows out of the steeper valley, around the mountains in the south-eastern corner of the domain. The modeled flow pattern stays nearly constant through the first eight hours of the s imulat ion, consistently over-estimating the easterly component of the flow in the eastern hal f of the domain, (see Figures 5.5.1-8a) This can be seen in the wind vector residual plots for these hours.(see Figures 5.5.1-8c) F igure 5.5.25 shows the values of the indices of agreement through the day of A u g u s t 23,1985, of two model runs , one wi th spatial ly var ied surface characteristics over the land and the other wi th spatial ly constant surface characteristics. The index of agreement of the run- made wi th spat ial ly varied surface characteristics starts out fa i r ly high at a value near .70 in the first few hours of simulation but then starts to decrease rapidly. The - 1 1 6 -decreasing agreement can be seen by examining the modeled and observed windfields for the hours between 04:00 and 08:00 when the observed flow is light and varied over most of the domain, (see Figures 5.5.4-8a and Figures 5.5.4-8b) The observed light and variable conditions are associated with the switch-over from the land to sea breeze regime. During this period the modeled flow is smooth over the domain and constant in time. It is not until 09:00 that the model starts to produce the switch. This gives poor agreement during the hours from 04:00 to 09:00.(see Figure 5.5.25) The main difference between the modeled and the observed would appear to be that the observed switch-over is preceded by a period of light and variable conditions where the model did not reproduce these same preceding conditions. As the sea breeze fills in, the agreement improves to a value between .6 and .7 for the hours between 10:00 and 19:00.(see Figure 5.5.25) During this period, the observed flow over the entire domain has a large westerly component to it. (see Figures 5.5.10-19b) The modeled flow on the other hand, stagnates over the eastern one third of the domain, (see Figures 5.5.10-19a) This westerly component of flow, not modeled in the eastern one third of the domain, may be caused by a large scale surface pressure gradient acting over the region containing the domain of the model. Intense surface heating in the interior of the province can create thermal lows that have effects reaching as far as the lower Fraser Valley. These surface lows also have much less diurnal variability than local pressure gradients created within the lower Fraser Valley. The evidence for this large westerly component of the flow being the result of a large scale pressure gradient becomes more convincing when the westerly component of the observed flow opposing the 850mb synoptic gradient, continues well into the night, (see Figures 5.5.19-22b) This is well after the time when one would expect the flow to have switched to a land breeze. Finally, throughout the entire day the wind vector plots show an over-prediction of the westerly component of the flow or conversely, an under-prediction of the easterly - 1 1 7 -component of the flow over the domain, (see Figures 5.5.1-24c) This is indicative of an effect acting over the whole domain. This adds to the above evidence for the existence of a large scale surface pressure gradient acting over the region. The model starts to produce the switch from sea to land breeze as early as 19:00 when the modeled surface heating goes to zero and the synoptic gradient begins to swing the flow back toward the west, (see Figure 5.5.19a) As the modeled land breeze fills in and the flow gains an easterly component over the entire domain, the agreement with the observed windfield decreases to its lowest value of 0.1 at midnight. This disparity between modeled and observed windfields could in part be caused by the above mentioned large scale surface pressure gradient. Another possible reason for this disagreement lies in the fact that modeled heat storage in the ground is very crudely parameterized and may not be representing the reality of the situation. This continued westerly component to the observed flow may, in part, be caused or enhanced by the positive fluxes of sensible heat created when the storage term in the energy budget (Qs) becomes negative. The release of stored heat from the earth surface may be available to counteract effects of radiative cooling at the surface. The release of stored heat is not modeled and therefore this effect, if it exists, is not modeled. The probable cause for the poor performance of the model on this day is a combination of the above factors. The thermal low in the interior of the province may be drawing air into the interior via the Fraser Valley. At the same time, the effects of radiative cooling may be balanced by storage losses. However, the latter is the smaller of the two effects and is of a shorter duration than the former. -118-Figures 5.5.1-24a These figures show the modeled windfields over the domain for the 24 modeled hours of August 23.19S5. The end of the one hour averaging periods are shown in the top left corner of the frames. The wind speed scales appear in the top right corner of the plot frames. Figures 5.5.1-24b These figures show the interpolated windfields for the 24 modeled hours of August 23,1985. The windfields shown are inverse distance weighted least squares fits of 2nd order polynomials to the data points. Only the sections of the windfields that are close to data collection sites are shown. This is because the reliability of the interpolation decreases rapidly with distance from measured points. Figures 5.5.1-24c These figures show the windfield residual plots for the 24 modeled hours of August 23,1985. The wind velocity residuals are the result of subtracting the observed from the predicted wind vectors. Circles enclosing the site numbers are plotted at the tip of the residual vectors with the tails of the residual vectors at the centres of the concentric circles. -119-RUG_23_I98S TIME (HRittTNl li 0 SURFACE WINDS tOOOfl CONTOUR 0.0 111.0 SCflLE (M/S) i 4) SO 80 DISTANCE (KM) — I — 90 100 Figure 5.5.1a 1985.08.23: 100 g U,,u I 1 I I 1 i — I * 0I3TPN?E IKMl Figure 5.5.1b - 1 2 0 -TIME. HOUR ENDING 1:00 RMSE (i/s) - 0.8925 RMSES («/s) = 0.5794 RMSEU («/s) - 0.6789 000 deg. 270 deg. D2 = 0.6445 HERN OBS. d/s) - 1.052 MEAN PRED. (n/s) - 0.88 SIGMA OBS. - 0.5769 SIGMA PRED. = 0.8721 090 deg. o.o s.o 5.0 SCALE (M/S) 180 deg. UIND VELOCITY RESIDUALS F i g u r e 5 . 5 . 1 c TIME, HOUR ENDING 2:00 RMSE (i/s) - 0.8348 RMSES («/s) = 0.4280 RMSEU («/s) - 0.7167 000 deg. 270 deg. 02 = 0.6930 MEAN OBS. (i/s) - 1.149 MEAN PRED. (n/s) - 0.91 SIGMA OBS. - 0.6192 SIGMA PRED. = 0.8954 090 deg. 0.0 3.0 6.0 SCALE (M/S) 180 deg. UIND VELOCITY RESIDUALS F i g u r e 5 . 5 . 2 c -121-Figure 5.5.2a 1985.08.23: 200 • i i u a 90 « 9_ • OISTONCE (KMI Figure 5.5.2b -122-Figure 5 . 5 . 3 b -123-TlflE. HOUR ENDING 3:00 RMSE In/s) - 0.9582 RMSES (i/s) = 0.5379 RMSEU («/s) - 0.7929 000 deg. 270 deg. 02 = 0.6120 HERN OBS. («/s) - 1.119 MEAN PRED. («/s) - 0.93 SIGMA OBS. - 0.6188 SIGMA PRED. = 0.9309 090 deg. 0.0 3.0 SCALE (M/S) 180 deg. i 8.0 UIND VELOCITY RESIDUALS Figure 5.5.3c TIME. HOUR ENDING 4:00 RMSE (i/s) - 1.1205 RMSES (i/s) = 0.6550 RMSEU («/s) - 0.9092 000 deg. 270 deg. 02 = 0.4101 MEAN OBS. (»/s) - 0.945 MEAN PRED. («/s) - 0.94 SIGMA OBS. - 0.4206 SIGMA PRED. = 0.9617 090 deg. — 0.0 3.0 6.0 SCALE (M/S) 180 deg. UIND VELOCITY RESIDUALS Figure 5.5.4c -124-Figure 5.5.4a -125-HJG_23_1985 TIrE (HR.WN) 5. 0 SURFACE WINDS iooon CONTOUR o.o 10.a SCPLE (M/S) 41 SO 80 OISTfiNCE (KM) Figure 5.5.5a 1985.08.23: 500 iooon CONTOUR 1 (M/?j * SfltfLED POINTS " * OISTPIZE (KM1 " Figure 5.5.5b - 1 2 6 -TIME, HOUR ENDING 5:00 RUSE d/s) - 1.1507 RMSES (»/s) =• 0.7808 RMSEU («/s) - 0.8453 000 deg. 270 deg. 02 = 0.3909 MEAN OBS. («/s) - 1.003 MEAN PRED. (»/s) - 0.91 SIGMA OBS. - 0.4549 SIGMA PRED. = 0.8803 090 deg. o.o 3.0 SCALE (M/S) 180 deg. UIND VELOCITY RESIDUALS F i g u r e 5 . 5 . 5 c TIME. HOUR ENOING 6:00 RMSE d/s) - 1.2003 RMSES («/s) = 0.7575 RMSEU In/s) - 0.9311 000 deg. 270 deg. 02 = 0.2904 MEAN OBS. («/s) - 1.091 MEAN PRED. d/s) - 0.92 SIGMA OBS. - 0.5586 SIGMA PRED. =• 0.9436 090 deg. 0.0 3.0 8.0 SCALE (M/S) 180 deg. UIND VELOCITY RESIDUALS F i g u r e 5 . 5 . 6 c - 1 2 7 -Figure 5.5.6b - 1 2 8 -F i g u r e 5 . 5 . 7 a - 1 2 9 -TIME. HOUR ENDING 7:00 RASE (i/s) - 1.2011 RMSES («/s) = 0.7733 RMSEU («/s) - 0.9191 270 deg. •2 « 0.2200 MERN OBS. d/s) - 1.034 MERN PRED. (»/s) - 0.89 SIGMR OBS. - 0.4960 SIGMR PRED. = 0.9357 000 deg. 090 deg. 180 deg. r SCRLE3'(M/S) W I N D V E L O C I T Y R E S I D U A L S F i g u r e 5 . 5 . 7 c TIME. HOUR ENDING 8:00 RMSE («/s) - 1.4185 RMSES <«/s) = 1.0791 RMSEU («/s) - 0.9208 000 deg. 270 deg. D2 = 0.4501 MERN OBS. (./si - 1.288 MERN PRED. In/s) - 0.89 SIGMR OBS. - 0.8575 SIGMR PRED. = 0.9357 090 deg. o.o s.o SCRLE (M/S) 180 deg. W I N D V E L O C I T Y R E S I D U A L S F i g u r e 5 . 5 . 8 c - 1 3 0 -F i g u r e 5 . 5 . 8 a 1985.08.23: 800 i i i F i g u r e 5 . 5 . 8 b - 1 3 1 -1985.08.23: 900 IU I 1 1 . - — ' • ' l . F i g u r e 5.5.9b -132-TlrlE. HOUR ENDING 9:00 RMSE d/s) - 1.3468 RMSES iu/s) = 1.2377 RMSEU («/s) - 0.5309 000 deg. 270 deg. 02 = 0.3664 MERN OBS. (i/s) - 0.864 MEAN PRED. d/s) - 0.46 SIGMA OBS. - 0.9304 SIGMA PRED. - 0.5495 090 deg. 0.0 3.0 s.o SCALE (M/S) 180 deg. UIND VELOCITY RESIDUALS Figure 5.5.9c TIME. HOUR ENDING 10:00 RMSE («/s) - 1.5711 RMSES (n/s) = 1.4900 RMSEU (»/s) - 0.4981 270 deg. D2 = 0.5122 MEAN OBS. («/s) - 1.460 MEAN PRED. In/s) - 0.45 SIGMA OBS. - 1.3638 SIGMA PRED. = 0.5505 000 deg. 090 deg. T 0.0 3.0 6.0 SCALE (M/S) 180 deg. UIND VELOCITY RESIDUALS Figure 5.5.10c - 1 3 3 -1985.08.23:1000 F i g u r e 5 . 5 . 1 0 b - 1 3 4 -1985.08.23:1100 Figure 5.5.11b - 1 3 5 -TIME. HOUR ENDING 11:00 RflSE d/s) - 1.5449 RMSES («/s) = 1.3720 RMSEU («/s) - 0.7102 000 deg. 270 deg. 02 = 0.5954 MERN OBS. («/s) - 1.675 MEAN PRED. («/s) - 0.72 SIGMR OBS. - 1.6175 SIGMR PRED. =• Q.8189 090 deg. o.o 3.0 SCALE (M/S) 180 deg. 6.0 WIND VELOCITY RESIDUALS Figure 5.5.11c TIME. HOUR ENDING 12:00 RMSE ti/s) - 1.3S22 RMSES (»/s) = 1.1333 RMSEU («/s) - 0.7375 270 deg. 02 = 0.6806 MEAN OBS. («/s) - 1.707 MEAN PRED. («/s) - 0.89 SIGMA OBS. - 1.6477 SIGMA PRED. = 0.9208 000 deg. 090 deg. 0.0 3.0 6.0 SCALE (M/S) 180 deg. WIND VELOCITY RESIDUALS Figure 5.5.12c -136-Figure 5.5.12a 1985.08.23:1200 Figure 5.5.12b - 1 3 7 -fiUG_23_1985 TIME (HRiMIN) 13J 0 SURFACE UINDS iooon CONTOUR I 1 o.g lo.i SCALE (M/S! 40 SO BO DISTANCE (KM) 70 00 90 100 Figure 5.5.13a 1985.08.23:1300 Figure 5.5.13b - 1 3 8 -TIME. HOUR ENOING 13:00 RMSE («/s) - 1.2333 RMSES l*/s) = 1.0115 RMSEU («/s) - 0.7056 000 deg. 270 deg. 02 = 0.7011 MEAN OBS. («/s) - 1.732 MEAN PREO. d/s) - 0.95 SIGMA OBS. - 1.5148 SIGMA PRED. = 0.8970 090 deg. ~i 1 0.0 3.0 s.o SCALE (M/S) 180 deg. UIND VELOCITY RESIDUALS F i g u r e 5 . 5 . 1 3 c TIME. HOUR ENDING 14:00 RMSE («/s) - 1.4263 RMSES («/s) = 1.2929 RMSEU («/s) - 0.6023 000 deg. 270 deg. D2 = 0.7108 MEAN OBS. («/s) - 2.024 MEAN PRED. (»/s) - 0.94 SIGMA OBS. - 1.7953 SIGMA PRED. = 0.8943 090 deg. 0.0 3.0 SCALE (M/S) 180 deg. UIND VELOCITY RESIDUALS F i g u r e 5 . 5 . 1 4 c -139-AUG_23J985 TIME (HRiHINl U i 0 SURFACE UINDS lOOOfl CONTOUR 0.0 10.0 SCALE (M/S) ' A U fa 40 SO 80 DISTANCE (KM) Figure 5.5.14a 1985.08.23:1400 i l l I ' , . - j r - l - ^ -Figure 5.5.14b -140-Figure 5.5.15a 1985.08.23:1500 Figure 5.5.15b - 1 4 1 -TIME. HOUR ENDING 15:00 RMSE («/s) - 1.8322 RMSES («/s) = 1.6977 RMSEU <«/s) - 0.6891 000 deg. 270 deg. D2 = 0.6282 MEAN OBS. («/s) - 2.358 MEAN PRED. («/s) - 0.92 SIGMA OBS. - 2.H50 SIGMA PRED. = 0.9081 090 deg. 0.0 3.0 SCALE (M/S) 180 deg. WIND VELOCITY RESIDUALS F i g u r e 5 . 5 . 1 5 c TIME, HOUR ENOING 16:00 RMSE d/s) - 1.8507 RMSES («/s) = 1.6891 RMSEU (»/s) - 0.7563 000 deg. 270 deg. 02 = 0.5973 MEAN OBS. (i/s) - 2.278 MEAN PRED. («/s) - 0.89 SIGMA OBS. - 2.0849 SIGMA PRED. = 0.9109 090 deg. - i 1 o.o 3.0 s.o SCALE (M/S) 180 deg. WIND VELOCITY RESIDUALS F i g u r e 5 . 5 . 1 6 c -142-Figure 5.5.16a 1985.08.23:1600 Figure 5.5.16b - 1 4 3 -1985.08.23:1700 i i t F i g u r e 5 . 5 . 1 7 b -144-000 deg. TIME. HOUR ENDING 17:00 RMSE (•/s) - 1.6498 090 deg. 180 deg. WIND VELOCITY RESIDUALS Figure 5.5.17c TIME. HOUR ENDING 18:00 RMSE («/s) - 1.4749 RMSES (n/s) = 1.3612 RMSEU ln/s) - 0.5677 000 deg. 270 deg. D2 = 0.6840 MERN OBS. («/s) - 1.798 MERN PRED. d/s) - 0.73 SIGMR OBS. - 1.5601 SIGMR PRED. = 0.8373 0.0 3.0 S.O SCALE (M/S) 180 deg. 090 deg. WIND VELOCITY RESIDUALS Figure 5.5.18c -145-Figure 5.5.18a 1985.08.23:1800 • ' ' „ Figure 5.5.18b -146-RUG_23_I985 H r E (HR.rUN) 1S> 0 SURFACE UINDS tooon CONTOUR o.o lo.o SCALE (M/S) 40 50 60 DISTANCE (KM) Figure 5.5.19a Figure 5.5.19b - 1 4 7 -TIME. HOUR ENOING 19:00 RMSE ( i / s ) - 113948 RMSES («/s) = 1.2889 RMSEU d /s) - 0.5331 000 deg. 270 deg. 02 = 0.6194 MEAN OBS. d /s ) - 1.431 MERN PRED. d /s) - 0.59 SIGMA OBS. - 1.2941 SIGMA PRED. = 0.6834 090 deg. T 0.0 S.O SCALE (M /S) 180 deg. l 6.0 UIND VELOCITY RESIDUALS F i g u r e 5 . 5 . 1 9 c TIME. HOUR ENDING 20:00 RMSE d /s) - 1.2640 RMSES ( i / s ) = 1.2107 RMSEU d /s ) - 0.3631 270 deg. 02 = 0.5532 MEAN OBS. d /s) - 0.874 MEAN PRED. ( i / s ) - 0.50 SIGMA OBS. - 1.0310 SIGMA PRED. " 0.4275 000 deg. 090 deg. - r 0.0 3.0 SCALE (M /S) 180 deg. I 6.0 UIND VELOCITY RESIDUALS F i g u r e 5 . 5 . 2 0 c - 1 4 8 -Figure 5.5.20a 1985.08.23:2000 Figure 5.5.20b -149-fiUG_23_198S TIME (HRjfUN) 21: 0 SURFACE UINDS loan CONTOUR 0.0 10.0 SCALE (M/S) 49 50 00 OISTRNCE (KM) Figure 5.5.21a 1985.08.23:210.0 Figure 5.5.21b - 1 5 0 -TIME. HOUR ENDING 21:00 RMSE (i/s) - 1.2916 RMSES d/s) = 1.2606 RMSEU d/s) - 0.2814 000 deg. 270 d e g . D2 = 0.4094 MERN OBS. d/s) - 0.625 MEAN PRED. d/s) - 0.72 SIGMA OBS. - 0.8367 SIGMA PRED. = 0.2853 090 deg. -i 1 0.0 3.0 6.0 SCALE (M/S) 180 deg. UIND VELOCITY RESIDUALS F i g u r e 5 . 5 . 2 1 c TIME. HOUR ENDING 22:00 RMSE d/s) - 1.2755 RMSES d/s) = 1.2008 RMSEU d/s) - 0.4301 270 d e g . D2 = 0.3366 MEAN OBS. d/s) - 0.447 MEAN PRED. d/s) - 0.82 SIGMA 08S. - 0.5710 SIGMA PRED. • 0.4541 000 deg. 090 deg. 0.0 3.0 6.0 SCALE (M/S) 180 deg. UIND VELOCITY RESIDUALS F i g u r e 5 . 5 . 2 2 c -151-AUG_23J985 TIME (HRiMIN) 22. 0 SURFACE WINDS lOOOfl CONTOUR o.o ip.o SCALE (M/S) — i r i 49 so ao DISTANCE (KM) Figure 5.5.22a 1985.08.23:2200 i if i ' , •. .'—. •„ Figure 5.5.22b -152-AUG_23_I985 TIME (HRtnTNJ 23. 0 SURFACE UINDS 1000M CONTOUR SCALE (M/s' a .a I <o so sa OISTRNCE (KM) 100 Figure 5.5.23a 1985.08.23:2300 iooon CONTOUR I 1 &a£ (M/S) * SWIED POINTS 01 5T»& (KM) Figure 5.5.23b - 1 5 3 -TIHE. HOUR ENOING 23:00 RMSE (i/s) - 1.4492 RMSES (n/s) = 1.3460 RMSEU <»/s) - 0.5371 000 deg. 270 deg. D2 = 0.1549 MERN OBS. («/s) - 0.520 MERN PRED. Im/s) - 0.86 SIGMA OBS. - 0.7105 SIGMA PRED. = 0.6344 090 deg. 0.0 3.0 S.O SCALE (M/S) 180 deg. UIND VELOCITY RESIDUALS F i g u r e 5 . 5 . 2 3 c TIME. HOUR ENDING 24:00 RMSE (»/s) - 1.5453 RMSES (»/s) =• 1.3812 RMSEU (n/s) - 0.6929 270 deg. D2 = 0.1207 MERN OBS. («/s) - 0.599 MERN PRED. l«/s) - 0.88 SIGMA OBS. - 0.7697 SIGMA PRED. = 0.7720 000 deg. 090 deg. - i 1 0.0 3.0 6.0 SCALE (M/S) 180 deg, UIND VELOCITY RESIDUALS F i g u r e 5 . 5 . 2 4 c -154-F i g u r e 5 . 5 . 2 4 a 1985.08.23:2400 i • F i g u r e 5 . 5 . 2 4 b INDICES OF AGREEMENT-SPRTIRLLY VARIED and SPHTIRLLY CONSTANT SURFACE CHARACTERISTICS (AUG.23/1985) 1- SPHTIRLLY VARIED SURFACE CHARACTERISTICS 2 - SPHTIRLLY CONSTANT SURFACE CHARACTERISTICS H—I—I-.-1—1—1-12 14 8 10 PACIFIC STANDARD TIME H—I—I—i 16 18 28 22 24 Ln I I 1 h—I 1 f—} ~ f -Figure 5.5.25 Indices of Agreement vs Time - August 23,1985 -156-5.6 G e n e r a l O b s e r v a t i o n s It seems reasonable to divide the lack of agreement into two distinct categories based on the cause of the disagreement. The first of the two is the frequent inability of the model to reproduce the timing of the switch from the land to the sea breeze and vice versa. It is during this transition period that the heating parameterization must be very accurate in order to force the model at the proper rate. With the intent of making the parameterization more accurate, specification of the spatial pattern of surface characteristics might be improved. This could be done by taking measurements of fluxes throughout the day at different locations in the domain. This would obtain a better estimate of the impacts of the spatial distribution of surface characteristics. A more simple and ultimately more realistic method would be experimentation with a number of patterns of surface characteristics, under a number of synoptic situations, in order to refine the patterns and magnitudes of the surface heating parameters. This might remove some of the systematic bias in the model and enable the model to reproduce the transitions in a more accurate fashion. The period during which the radiative balance forces values of surface heating from negative to positive and vice versa, seems to be the most critical from a modeling point of view. It is during this time that the height of the thermally affected boundary layer is changing the most. The inability of the model to account for these changes is certainly one of the factors affecting the accuracy of the model during these periods (Mass and Dempsey, 1985b). This could in part be remedied by the parameterization of the spatial and temporal variability of the thermally affected boundary layer. This could be done without drastically changing the formulation of the model if a physically realistic parameterization scheme were derived. The second and the more fundamental reason for the seemingly poor performance of the model is embodied in the method by which it is validated. It is tempting to look at the plots of the index of agreement through the day and assess the worth of the model to be - 1 5 7 -very low. This would be a premature assessment without further consideration of what the model is up against in terms of the task it is being required to perform and the resolution at which it operates. The measurements used to validate the model are point measurements taken at i rregular ly spaced intervals over a l imited section of the domain. These measurements are taken in flow that is affected by objects that range in size from as smal l as an individual roughness element to as large as the mountains bordering the domain (Taylor and Lee, 1984; H u n t and Simpson, 1982) Embedded w i t h i n the flow are eddies that are mechanically and thermal ly forced and that v a r y tremendously in scales of length and time. A great deal of this variabi l i ty is eliminated by time averaging but this sti l l leaves substantial spat ia l var iabi l i ty to account for. The model has a resolution of roughly four kilometres and a time step of 60 seconds. A n y t h i n g that happens in the real atmosphere at length or time scales shorter than these w i l l not be represented except in a very approximate w a y . It seems unreasonable to expect the model to produce agreement to this scale of resolution when it is real ly only designed to simulate a very general flow pattern. A s mentioned earlier in Section 4.2, the model cannot be expected to work at Froude numbers much greater than 1.0. This sensit ivity to Froude number is a l imitat ion that can be used to demonstrate how the model can produce poor statistical agreement while doing a good job of modeling the general flow pattern. Calculations presented in Section 4.2 show that U / h , the wind speed divided by the height of the barr ier , should not exceed 10 s for the modeled flow to be properly deflected around an obstacle. It can be seen that as the height of the obstruction decreases, this l imit is approached. The flow around obstacles wi th a height of up to a hundred metres would possibly not be modeled even if the obstacles were horizontally larger than the four kilometre resolution of the model. This brings to mind the question of the degree to which the model actual ly "feels" the topography. -158-Figures 5.6.1 and 5.6.2 are plots of the coefficient of variation against the mean wind speed for each of the 24 hours modeled. The coefficient of variation is the standard deviation of the wind speed over the domain divided by the mean wind speed over the domain. One would expect, given the model's lack of ability to "feel" the topography at higher Froude numbers and the fact that the Froude number increases with wind speed, the spatial variation and therefore the coefficient of variation would decrease as wind speed increases. This trend is clearly demonstrated by the observed values on July 20 and, although not as clearly, can also be seen on August 23. In the modeled values this trend is only slightly visible to the optimistic eye on August 23 and does not appear at all on July 20. This is an indication that the model is operating near the limit of where it can be expected to perform. This limit is imposed by the Froude number of the situations being modeled. This is of fundamental significance in terms of what can be expected of the model in an objective validation using the observation network employed in this study. All of the measurement sites were on the valley floor where a great percentage of the perturbations in the flow caused by small variations in terrain height, could be measured but not modeled. Validation sites nearer the mountains might have given the model more credit for producing channeling that is clearly visible in the modeled windfields. This would also give more favourable statistical values in the validation process. z O 1 n 1-<r i—i > Lu ° . 6 -h-Z u n U. U. O • 2 <J 1 n COEFFICIENT OF VflRinTION v s NIND SPEED Q JUL.20/1985 U © OBSERVED QQ. Q ^ MODELED 0 © 0 © © 0 0 0 00 Q 2< IS IS 1' it3 22 2a a 21 1 1 1 1 1 1 1 I I -) 0 t —t " "I" ~- 1 1 I I 1 1 1 1 I 2 3 4 WIND S P E E D ( m / s ) Figure 5.6.1 Coefficient of Variation vs Wind Speed - Ju l y 20,1985 c 1-1 fB. U l O N COEFFICIENT OF VRRIRTION o o fD i-h Hi H-O H-fD 3 3 O o cn a in n n I- 70 m < • n • n o m ~n i—i n t—t n < i-i to rt H* O 3 < 3 CL CO -a ft> (B a. > c C cn 00 U l cn TJ m a 3 \ o < D 73 ID JJ C L7) ru LO \ CD CD (_n CD TJ m m a - 0 9 1 -- 1 6 1 -6. CONCLUSIONS AND RECOMMENDATIONS FOR FURTHER WORK The objectives of this study were to examine the ability of a one-level mesoscale numerical model to simulate thermally and mechanically forced flows over complex terrain. This attempt met with mixed success in that the model did produce a number of the more salient features of these flows and circulations, though the finer details are in some cases not well modeled. The model, as used in this study, responds well to spatial variations of surface heating and cooling but is somewhat less sensitive to the spatial variability of surface roughness. It was found that changes in the specification of the spatial pattern of the surface heating and cooling created substantial changes in modeled windfields throughout the day. These changes had a substantial affect on the agreement between modeled and measured windfields. When used in situations where the domain is under the influence of a weak synoptic gradient, the local thermal forcing was found to dominate. This implies that the specification of the upper boundary conditions is relatively unimportant in the presence of a slack synoptic gradient. Small scale variations in the observed flow caused by small scale topographic and roughness features were not well represented by the model. Part of the inability of the model to reproduce these small scale variations originates in the model formulation and is caused by attempting to use the model in situations where the Froude number is too high. Calculations showed that the model could only be expected to reproduce the larger scale channeling effects while leaving the smaller scale variations in the flow observed over the valley bottom, unrepresented in the modeled windfields. In the cases presented here, use of the model at high Froude numbers was done at the cost of agreement between measured and modeled windfields. - 1 6 2 -Use of the model is clearly limited in certain situations including those where complex vertical thermal structure exists or in high Froude number situations (Mass and Dempsey, 1985b). In modeling thermally driven, mechanically forced flows and circulations, it has been shown that the latter of the two limitations comes into play. This fundamental limitation ultimately places restrictions on the ability of the model to produce accurate, verifiable results at the resolution used in this study. There are a number of changes that might be made to the model in order to improve its performance. Mass and Dempsey (1985b) suggest that a number of changes might be made to the model to improve its performance. These include: i) a parameterization of the depth of the layer of topographic influence; ii) parameterizing the rising and sinking of air caused by low level convergence and divergence; iii) improving the parameterization of friction to make it stability dependent and spatially variable; iv) allowing the upper atmosphere lapse rate to change spatially; and v) including a better parameterization of diabatic effects by making the depth to which the atmosphere is affected by surface heating, variable in time and space. In using the one-level model made available by Mass and Dempsey, a number of changes, including some of the above, were implemented. Specifically an attempt was made to include spatial variability in the diabatic forcing and friction parameterizations. As is implied in the text of this thesis, the phenomenon which is modeled is one of the major factors in determining model formulation. In modeling phenomena such as thermally induced flows over terrain similar to the lower Fraser Valley of British Columbia, there are a number of further directions that might be taken in order to improve the performance of the model. The first of these directions involves more the improvement of surface characteristic specification than model formulation itself. More accurate specification of spatial and temporal variation of surface heating in combination with a -163-sensitivity analysis under varied synoptic situations, could yield a significant improvement in performance. This refinement could be accomplished by either using measurement of fluxes over the real domain to determine the pattern of surface heating or by satellite-based assessment of spatial variability of surface fluxes based on the methods of Carlson et al. (1981). Another method might be the use of numerical experimentation to determine the pattern that yields the best agreement. Another factor that affects model performance within the context of this study is the depth to which the atmosphere is affected by surface heating. It is known that the depth of the mixed layer over the domain and therefore the sensitivity of the surface pressure to surface heat fluxes, varies spatially and temporally (Steyn and Oke, 1982). Therefore, as stated by Mass and Dempsey (1985b), improvment in model performance may be realized if the depth to which the atmosphere is affected by diabatic forcing at the surface, is made variable by parameterizing a mixed layer that is both spatially and temporally variable. The effects of large scale thermally forced surface pressure gradients have been shown to cause poor model performance in that they are not represented by the model. A method by which the agreement of the model could be enhanced is the inclusion of the present model domain within a much lower resolution grid. This would allow the modeling of larger scale forcing which has shown to cause disparity between modeled and observed fields. The resolution of the model might be increased in an attempt to increase agreement with the observed fields. However, this may not be a wise choice of strategies. With a resolution of approximately four kilometres over the domain, and using a three grid space buffer at the edges, the model uses some 850 seconds of computer time to do a 24 hour simulation. This is a ratio of approximately l:100-simulated to real time. A fully three-dimensional model, run with the same resolution over the same domain has a 1:4-simulated to real time ratio (McKendry, 1986). In comparing the two types of models, the one-level model when used in an appropriate situation is, by a great margin, more effective at producing a surface windfield prediction. If the resolution of the one-level model is - 1 6 4 -doubled, the time step must be reduced to roughly half the present value with the result that the model will take eight times as long to run on the same computer. This makes the use of the one-level model much less attractive when it is compared to the three-dimensional model which in most cases gives better results. At the same time, the increase in resolution may not improve the performance of the model if it is being used near the limiting Froude number of 1.0. The method of validation, although fundamentally sound, is not well matched to the model structure. The discretization of real fields within the modeling process imparts a certain amount of averaging over the fields. This spatial and temporal averaging is not well represented in the validation process. This may be a drawback that must be considered inherent within a validation procedure of this nature. To match the temporal and spatial averaging done by the model would require great amounts of equipment and effort. Even with the availability of the expertise involved in such an undertaking, the cost of the equipment required would certainly prohibit such an endeavor. The statistical and graphical methods used within the validation process of this study were found very useful in assessing model performance as well as within a diagnostic framework. These methods also make possible the objective comparison of this work with other modeling attempts. The one-level mesoscale model is a simple, economic tool for diagnosing surface winds over complex terrain. It has some limitations that hinder performance in certain situations. In producing the general flow pattern at the surface, the performance of the model is satisfactory under the conditions used in the simulations presented here. If a more detailed flow pattern is required in high Froude number situations, the use of this model, even with its improved parameterization, is inappropriate. -165-BIBILIOGRAPHY Anderson, G. E. (1971) Mesoscale Influences on Windfields Journal of Applied Meteorology, Vol. 10, p. 377-386 Anderson, D. A., J. C. Tannehill and R. H. Pletcher (1984) Computational Fluid Mechanics and Heat Transfer McGraw-Hill, New York, 599pp. Angell, J. K. and D. H. Pack (1965) Study of the Sea Breeze at Atlantic City, N. J. using Tetroons as Lagrangian Tracers Monthly Weather Review, Vol. 93, p. 475-493 Atkinson, B. W. (1981) Mesoscale Atmospheric Circulations Academic Press, London, England, 495pp. Atkinson, B. W. (1984) The Mesoscale Atmosphere Inaugral Lecture, Dept. of Geography and Earth Science, Queen Mary College, University of London Ayotte, K. W. (1985) Adjustment to One Elevation of Wind Speed Measurements made at several Elevation Unpublished Manuscript. Beer, T. (1974) Atmospheric Waves John Wiley and Sons, Toronto, 300pp. Bornstein, R. D. and D. S. Johnson (1977) Urban-Rural Wind Velocity Differences Atmospheric Environment, Vol. 11, p. 597-604 Burden, R.L., J.D. Faires and A.C. Reynolds (1981) Numerical Analysis. 2nd Edition Prindle, Weber and Schmidt, Boston, 598pp Carlson, T. N. et. al. (1981) Satellite Estimation of the Surface Energy Balance, Moisture Availability and Thermal Inertia Journal of Applied Meteorology, Vol. 20, No. 1, p. 67-87 - 1 6 6 -Cleugh, H. A. and T. R. Oke (1986) Suburban-Rural Energy Balance Comparisons in Summer for Vancouver, B. C. Boundary Layer Meteorology, In Press Concord Scientific Corporation (1982) Vancouver Oxidants Study - Air Quality Study - Final Report prepared for Environmental Protection Service, Environment Canada Concord Scientific Corporation (1985) Vancouver Oxidants Study - Air Quality Analysis Update, 1982-1984 - Final Report prepared for Environmental Protection Service, Environment Canada Counihan, J. (1975) Adiabatic Atmospheric Boundary Layers: A Review and Analysis of Data from the Period 1880-1972 Atmospheric Environment, Vol. 9, p. 871-905 Danard, M. B. (1971) A Numerical Study of the Effects of Longwave Radiation and Surface Friction on Cyclone Development Monthly Weather Review, Vol. 99, p. 831-839 Danard, M. B.. (1977) A Simple Model For Mesoscale Effects of Topography on Surface Winds Monthly Weather Review, Vol. 105, p. 572-581 Danard, M. and B. Thompson (1983) Modeling Winds in Lancaster Sound and Northwestern Baffin Bay Atmosphere-Ocean, Vol. 21, No. 1, p. 69-81 Deardorff, J. W. (1972) Parameterization of the Planetary Boundary Layer for use in General Circulation Models Monthly Weather Review, Vol. 100, p. 93-106 Dickerson, M. H. (1978) MASCON - A Mass Consistent Atmosphere Flux Model for Regions of Complex Terrain Journal of Applied Meteorology, Vol. 17, p. 241-253 Endlich, R. M. (1967) An Iterative Method of Altering the Kinematic Properties of Wind Fields Journal of Applied Meteorology, Vol. 6, p. 837-844 -167-Estoque, M. A. (1973) Numerical Modeling of the Planetary Boundary Layer in Workshop on Micrometeorology, American Meteorological Society Publication, Boston, Mass., p. 217-270 Findlater, J. (1963) Some Aerial Explorations of Coastal Airflow Meteorological Magazine, Vol. 92, p. 231-243 Fisher, E. L. (1960) An Observational Study of the Sea Breeze Journal of Meteorology, Vol. 17, p. 645-660 Frizzola, J. A. and E. L. Fisher (1963) A Series of Sea Breeze Observations in the New York City Area Journal of Applied Meteorology, Vol. 2, p. 722-739 Gerrity, J.P., R.D. McPherson and P.D. Polger (1972) On the Efficient Reduction of Truncation Error in Numerical Prediction Models Monthly Weather Review, Vol. 100, p. 637-643 Goodin, W. R., G. J. McRae and J. H. Seinfeld (1979) A Comparison of Interpolation Methods for Sparse Data: Applications to Wind and Concentration Fields Journal of Applied Meteorology, Vol. 2,.p. 761-771 Haltiner, G. J. and R. T. Williams (1980) Numerical Prediction and Dynamic Meteorology John Wiley and Sons Inc., Toronto, 477pp. Haurwitz, B. (1947) Comments on the Sea Breeze Circulation Journal of Meteorology, Vol. 4, p. 1-8 Hay, J. E. and D. C. McKay. (1985) Estimating Solar Irradiance on Inclined Surfaces: A Review and Assessment of Methodologies International Journal of Solar Energy, Vol. 3, p. 203-240 Hay, J. E. and T. R. Oke. (1976) The Climate of Vancouver B. C. Geographical Series, No. 23, Tantalus Research Ltd., Vancouver, 48pp. - 1 6 8 -Hjelmfelt, M. R. (1982) Numerical Simulation and the Effects of St. Louis on Mesoscale Boundary-Layer Airflow and Vertical Air Motion: Simulations of Urban vs Non-Urban Effects Journal of Applied Meteorology, Vol. 21. No. 9, p. 1239-1257 Holton, J. R. (1979) An Introduction to Dynamic Meteorology Academic Press Inc., New York, 391pp. Hunt, J. C. R. and J. E. Simpson (1982) Atmospheric Boundary Layers over Non-homogeneous Terrain in Engineering Meteorology. E. Plate, ed. Elsiver Publishing, New York, p. 269-314 Idso, S. B. and R. D. Jackson. (1969) Thermal Radiation from the Atmosphere Journal of Geophysical Research, Vol. 74, p. 5397-5403 Kalanda, B. D. (1979) Suburban Evaporation Estimates in Vancouver, B. C. M.Sc. Thesis. University of British Columbia, Vancouver Keyser, D. and R. A. Anthes. (1977) The Applicability of a Mixed-Layer Model of the Planetary Boundary Layer to Real-Data Forecasting Monthly Weather Review, Vol. 105, p. 1351-1371 Kimble, G. H. T. et. al. (1946) Tropical Land and Sea Breezes Bulletin of the American Meteorological Society, Vol. 27, p. 99-113 Lavoie, R. L. (1972) A Mesoscale Model of Lake-effect Storms Journal of the Atmospheric Sciences, Vol. 29, p. 1025-1040 Lavoie, R. L. (1974) A Numerical Model of Trade Wind Weather on Oahu. Monthly Weather Review, Vol. 102, p. 630-637 Leblond, P. H. and L. A. Mysak (1978) Waves in the Ocean Elsevier Scientific Publishing Company, New York 602pp. - 1 6 9 -Mass. C. F. (1981) Topographically Forced Convergence in Western Washington State Monthly Weather Review, Vol. 109, p. 1335-1347 Mass, C. F. (1982) The Topographically Forced Diurnal Circulations of Western Washington State and their Influence on Precipitation Monthly Weather Review, Vol. 110, p. 170-183 Mass, C. F. and D. P. Dempsey. (1985a) A Topographically Forced Convergence Line in the Lee of the Olympic Mountains Monthly Weather Review, Vol. 113, p. 659-663 Mass, C. F. and D. P. Dempsey. (1985b) A One-Level, Mesoscale Model for Diagnosing Surface Winds in Mountainous and Coastal Regions Monthly Weather Review, Vol. 113, p. 1211-1227 McKendry, I. G. (1986) Personal Communication McKendry, I. G. and D. G. Steyn (1986) Application of a 3-D Mesoscale Model to Windfield Prediction in the Lower Fraser Valley, B. C. Presented to the 20th Annual Congress of the Canadian Meteorologic and Oceanographic Society, Regina. McPherson, R. D. (1970) A Numerical Study of the Effect of a Coastal Irregularity on the Sea Breeze Journal of Applied Meteorology, Vol. 9, p. 767-777 Mesinger, F. and A. Arakawa (1976) Numerical Methods used in Atmospheric Models GARP Publishing Series No. 17, Vol. 1, 64pp. Neumann, J. (1977) On the Rotation Rate of the Direction of Sea and Land Breezes Journal of Atmospheric Sciences, Vol. 34, p. 1913-1917 Oke, T. R. (1978) Boundary Layer Climates Methuen and Co. Ltd., London, 372pp. -170-Oke, T . R., B. D. Kalanda and D. G. Steyn. (1981) Parameterization of Heat Storage in Urban Areas Urban Ecology, Vol. 5 (1980/1981), p. 45-54 Overland, R. E. and B. A. Walter, Jr. (1981) Gap Winds in the Strait of Juan de Fuca Monthly Weather Review, Vol. 109, p. 2221-2233 Overland, J. E. and B. A. Walter, Jr. (1983) Marine Weather of the Inland Waters of Western Washington NOAA Technical Memorandum ERL PMEL-44. 62pp. Panofsky, H. A. and J. A. Dutton. (1984) Atmosphere Turbulence: Models and Methods for Engineering Applications John Wiley and Sons Inc., Toronto, 397pp. Pasquill, F. and F. B. Smith. (1983) Atmospheric Diffusion Ellis Horwood, New York, 437pp. Pielke, R. A. (1984) Mesoscale Meteorological Modeling Academic Press Inc., Orlando, Florida, 612pp. Sasaki, Y. (1970) Some Basic Formalisms in Numerical Variational Analysis Monthly Weather Review, Vol. 98, p. 875-883 Schoenberg, S. A. (1983) Regional Wind Patterns of the Inland Waters of Western Washington and Southern Britisn Columbia NOAA Technical Memorandum ERL PMEL-43 Sherman, C A. (1978) A Mass-Consistent Model for Wind Fields Over Complex Terrain Journal of Applied Meteorology, Vol. 17, p. 312-319 Silversides, R. H. (1978) Forest and Airport Wind Speeds Atmosphere-Ocean, Vol. 16, No. 3, p. 293-299 - 1 7 1 -Smagorinsky, J. (1981) Epilogue: Perspective of Dynamical Meteorology in Dynamical Meteorology. An Introductory Selection. B.W. Atkinson, ed. Methuen, p. 205-219 Spencer, J . W. (1971) Fourier Series Representation of the Position of the Sun Search, Vol. 2, p. 172 Steyn, D. G. (1976) Computation of Azimuths, Slope Angles and Surface Normals over a given Topograph}' The South African Geographical Journal, Vol. 58, No. 2, p. 130-133 Steyn, D. G. (1980) Turbulence. Diffusion and the Daytime Mixed Layer over a Coastal City Ph. D. Thesis. University of British Columbia, Vancouver Steyn, D, G. (19S4) Observational Studies of Sea Breezes in the Lower Fraser Valley, B. C : Preliminary Results Presented to the 18th Annual Congress of the Canadian Meteorologic and Oceanographic Society, Halifax Steyn, D. G. and D. A. Faulkner, D. A. (1986) The Climatology of Sea Breezes in the Lower Fraser Valley, B. C. CMOS Climatological Bulletin, Vol. 20, No. 3, in press Steyn, D. G. and I. G. McKendry (1986) Sea Breeze Dynamics in the Lower Fraser Valley, B. C. Research Paper under contract with B.C.Ministry of Environment and the University of British Columbia, 41pp. Steyn, D. G. and T. R. Oke (1982) The Depth of the Daytime Mixed Layer at Two Coastal Sites: A Model and its Validation Boundary-Layer Meteorology, Vol. 24, p. 161-180 Sutcliffe, R. C. (1937) The Sea Breeze at Felixstowe; A Statistical Investigation of Pilot Balloon Ascents up to 5,500 feet Quarterly Journal of the Royal Meteorological Society, Vol. 63, p. 137-146 - 1 7 2 -Taylor, P. A. and R. J. Lee (1984) Simple Guidelines for Estimating Wind Speed Variations due to Small Scale Topographic Features CMOS Climatological Bulletin, Vol. 18, No.2, p. 3-32 Theissen, A. H. (1911) Precipitation Averages for Large Areas Monthly Weather Review, p. 1082-1084 Vasanji, Z. and I. S. Gartshore (1978) An Analysis of Six Simultaneous Wind Records Taken Near Vancouver Atmosphere-Ocean, Vol. 16, No. 2, p. 145-156 Wexler, R. (1946) Theory and Observations of Land and Sea Breezes Bulletin of the American Meteorological Society, Vol. 27, p. 272-287 Willmott, C. J. (1981) On the Validation of Models Physical Geography, Vol. 2, p. 184-194 Willmott, C. J. (1984) On the Evaluation of Model Performance in Physical Geography in Spatial Statistics and Models, G. L. Gaile and C. J. Willmott, eds., D. Reidel Publishing Co., Hingham, Mass., p. 443-460 Willmott C. J. et. al. (1985) Statistics for the Evaluation and Comparison of Models Journal of Geophysical Research, Vol. 90, No. C5, p. 8995-9005 Wilson, R. G., J. B. Mills and E. P. Wituschek (1984) A Report on the Assessment of Photochemical Oxidants in the Lower Mainland Queens Printers, Victoria, 23pp. Wippermann, F. (1981) The Applicability of Several Approximations in Mesoscale Modeling: A Linear Approach Contributions to Atmospheric Physics, Vol. 54, No. 2, p. 298-308 -173-A PPENDIX I: S1GMA-C00RD1NATE TRANSFORMATION The formulation of the model equations makes use of a e-coordinate transformation in order to simplify the governing equations over a domain which has varying surface elevation. This coordinate transformation operates only on the vertical coordinate while leaving the horizontal coordinates unaffected. The transformed vertical coordinate is: The reason for its usefulness is that it is a terrain following coordinate, ie. the value of Sigma is constant at all points on the earth's surface. This is particularly useful in the case of a one-level model because the only level at which the governing equations need to be formulated is the surface where a has the value of unity. The transformation affects the vertical variable and therefore any of the derivatives, with respect to the vertical coordinate. Stated more precisely: a = P/Ps eq.al. 1 where P is the pressure at elevation and P g is the surface pressure. A(x,y,o,t) = A(x,y,z(x,y,o,t),t) eq.al.2 If a partial derivative is taken with respect to s where s = x,y,t , the result is: eq.al.3 In turn, the transformed version of the grad and the grad-dot operators become: -174-V QA = V ZA + !A l £v az V a-A = V z-A + 1**1 V A Z eq.al.4 3o 3Z 3a 3z (Haltiner and Williams, 1980) The above transformation is used in both the velocitj' tendency equation and the temperature tendency equation. The example presented here is the result of the transformation operating on the grad-Pg terms in the model velocity tendency equation applied at the surface. 1^ 1 = -V s*v aV s - f k x V s - { g v a z s + R T s v a l n P s } + F + K M v a 2 v s eq.al.5 d t ...as presented bj' Mass and Dempsey(1985b) To illustrate the effect of the transform in a graphical fashion, Figures a l . l and al.2 show the two vertical coordinate systems. Figure a l . l shows surfaces of constant pressure above a hypothetical topography. For the sake of demonstration the atmosphere is an ideal one that decays in pressure logrithmically with Z. The dashed line represents the topography and the heavy solid line which varies sinusoidally across the domain is the datum pressure with the scale to the left. Surfaces of constant pressure are plotted with pressures ranging from lOOOmb to 850mb in 50mb increments. It is important to note that at the surface the pressure varies substantially. This makes the use of pressure coordinates a less than simple system to work with. Figure al.2 shows the same topograph}' and atmosphere. In this figure, surfaces of constant e are plotted. It can be seen that these surfaces follow the earth's surface. The « f = 1 surface and the earth's surface coincide. As mentioned previousl}', this is a tremendous advantage when working within a domain that has a varying lower surface elevation. 1803 F i g u r e a l . l Constant P r e s s u r e S u r f a c e s Over Topography - 9 Z T --177-APPENDIX II: DERIVATION OF THE TENDENCY EQUATIONS Temperature Tendency The derivation of the temperature tendency equation starts with a simple statement of the first law of thermodynamics: Cpil - c u L = Q eq.a2.1 at at where Q is heat input. Note that the time derivatives of temperature and pressure are material derivatives. The total derivatives are expanded to show the local and advective changes of temperature with time. Cpfll + o i l + vlL + v^ I} - «{E + uiL + v— + wil} = Q eq.a2.2 at ax ay az at ax ay 3z The vertical terms are separated from the horizontal terms: — The temperature tendency equation is then transformed into sigma-coordinates and evaluated at the surface where sigma is unity everywhere over the domain at all times. The actual transformation process is described in Appendix I. r 3 T c 3 T C 3 T c i f3Pc , 3Pc , a P o r 3 P, - i + u — i + v—§-) - C n l — 5 - + u—2. + v — i f - a { — i + V H * V 0 P S } at ax ay p at ax ay at eq.a2.4 { V H i f l V a z s } + { % V a z s C f i I s } - {V> V 0 z s a i f i } = Q 3 Z o Z o Z Cancellation of the grad-Z terms leaves; -178-cp{ill + V H VaTs} -a{^l + V R VaPs} = Q eq.a2.5 9 t 31 ^here : a = \/0 = pj/p Further manipulation and the addition of a diffusion term yeilds; 3T i = _vs V a T s + E s j U ^ S + v s V a l n P s } + i - + K TV H 2 T S e ^ 2 - 6 3t Cp 3t Cp ...as presented by Mass and Dempsey(19S5b). Note that the grad n-T term is calculated as a horizontal gradient. If this gradient is calculated on a surface holding sigma constant, then the variation in elevation over the surface where sigma equals 1.0 will induce errors in the diffusion term calculation. The horizontal diffusion term is calculated by adjusting the temperatures to one level using the modeled lapse rate and then calculating the horizontal gradient of those temperatures. Velocity Tendency Newton's second law applied to atmospheric motion is the following: dV _ = -aVP + q dt y eq.a2.7 As with the temperature tendency equation, the derivative involved is a total derivative. Expansion of the total derivative as well as the inclusion of terms to account for the effect of a rotating reference frame and friction gives the following; - 1 7 9 -I V + v . v V = „ a V p . 2 n x v + i + ? eq.a2.8 Applying the equation at the surface, with the inclusion of a diffusion term and the use of the coordinate transform described in Appendix I, equation a2.8 can be rewritten in sigma-coordinates: — = -V S«V 0V S - f k x V s - { g V a z s + RT sV alnP s} + F + KMV a 2V s eq.a2.9 3 t as presented by Mass and Dempsey(1985b) -180-LIST OF ANEMOMETER LOCATIONS Number Location Administered by 1 2294 West 10th, Kitsilano G.V.R.D. 2 6200 Block, East Hastings St. G.V.R.D. 3 Confederation Park G.V.R.D. 4 75 Riverside Drive, N. Van. G.V.R.D. 5 Anmore G.V.R.D. 6 Rocky Point Park G.V.R.D. 7 South Langley B.C.M.O.E. 8 Sumas B.C.M.O.E. 9 Surrey B.C.M.O.E. 10 Langley U.B.C. 11 Boundary Bay U.B.C. 12 Cloverdale U.B.C. 13 Aldergrove U.B.C. 14 Pitt Meadows U.B.C. 15 South Arm Park, Richmond U.B.C. 16 White Rock U.B.C. 17 Masterwash Products Ltd. U.B.C. 18 Abbottsford A.E.S. 19 Vancouver Airport A.E.S. 20 Jericho Yacht Club A.E.S. 21 Tsawwassen A.E.S. 22 Vancouver Harbour A.E.S. 23 Cranberry Farm U.B.C. 24 Sand Heads A.E.S. 25 Point Atkison A.E.S. Table a3.1 List of locations of instrument sites used in validation of model, (after Steyn and McKendry, 1986) o I 1 1 1 1 1 i i 1 i I 0.0 10.0 20.0 30.0 10.0 50.0 60.0 70.0 80 .0 90.0 100.0 (KM) Figure a3.1 Validation Data Collection Sites (after Steyn and McKendry,1986) -182-APPENDIX III: COLLECTION OF VALIDATION DATA In order to validate model results, measurements of wind speed and direction were made over the model domain. A list of these stations and their locations appears in Table a3.1. The locations are plotted on a map of the domain in Figure a3.1. These stations all belong to one of four groups: those operated by the Atmospheric Environment Service (A.E.S), those operated by the Greater Vancouver Regional District (G.V.R.D.), those operated by the British Columbia Ministry of Environment (B.C.M.O.E.) and those operated by a group doing research under the direction of Dr. D.Steyn at the Department of Geography at the University of British Columbia. With the.exception of the latter group, the quality control and collection of data was out of the hands of the author. The wind speed data collected at the sites operated by the G.V.R.D are somewhat suspect for reasons outlined by Ayotte (1985). The effects of building wake and the inhomogeneity and inadequate length of fetch have led the author to treat the data collected at these sites with a large degree of skepticism. For this reason, the data from these station are not included in the statistical treatment of the modeled fields. A. E.S. STATIONS The A.E.S. stations all conform to World Meteorological Organization (W.M.O.) standards for site specification of wind measurement and all of the measurements were taken at 10 metres above the local surface. The instruments used are of the U2A type. These instruments sample continuously and the continuous data were averaged hourly at all of the stations used. B. C.M.O.E. STATIONS The B.C.M.O.E. stations are all located within areas of adequate homogeneous fetch and, like the A.E.S. stations, all measurements were taken at 10 metres above the -183-surface. The anemometers used were Anderaa anemometers with Gill Microvane wind direction sensors. The anemometers had a stall speed of 0.5 m/s. The data were sampled every 10 seconds and stored on Campbell Scientific CR21 data loggers. The speed and direction data were vector averaged over one hour periods. U.B.C STATIONS The U.B.C. stations were located at ten sites over the lower Fraser Valley. The site selections were made in such a way as to complement the already existing A.E.S. and B.C.M.O.E. network. Because of concerns for security and the need for the cooperation of landowners, local, provincial and federal governments whose properties the instrumentation was to be installed upon, the site choice could not always be optimal. Fortunately, fetch requirements were met in all cases. The instruments used to measure wind speed and direction at all of the U.B.C. sites were Met-One wind speed and direction sensors. These were linked to Campbell Scientific CR21 data loggers in all cases. The sampling frequency was six times per minute and the averaging time was ten minutes. The wind speed and direction were broken into vector components and a vector average was performed. At the Pitt Meadows, Aldergrove C.F.S., Cloverdale and Cranberry Farm sites the instruments were mounted on a 1 metre long cross-arm at the top of a 10 meter length of aluminum tubing. These towers were supported by two sets of three stays radiating out from points at the top and two-thirds of the way up the tower. It was not always possible or desirable to use the towers described above so various other methods of mounting on a number of different towers were used. The instrumentation at the Langley and Boundary Bay Airports was affixed to 2 meter cross-arms at the top of already existing A.E.S. instrument towers, parallel to already existing instruments. The reason the A.E.S. instruments were not used was because the instrumentation was not of the recording type. By virtue of the fact that these - 1 8 4 -sites were at airports, they had large open areas surrounding them that were free of high obstacles. This setting provided excellent fetch in all directions. The Masterwash site instrumentation was mounted on a 2 metre cross-arm attached to the top of a radio antenna on top of the Masterwash Corporation headquarters in the Tilbury Industrial Park in the municipality of Delta. The building is three stories (10 metres) in height and would be cause for concern pertaining to wake effects on wind data had the radio tower not been slightly more than 10 metres high itself. This put the instruments well above the region of accelerated flow that is commonly associated with flow over bluff bodies. The South Arm Park site was at a neighborhood swimming pool in suburban Richmond. The instrumentation was affixed to a 2 metre cross-arm at the 10 metre level of a lamp post at the pool's edge. The surface surrounding the location was similar in roughness in all directions with greatest part of the collection of roughness elements being one and two storey dwellings in a suburban setting. The White Rock installation was on top of the White Rock City Hall. Again the instruments were mounted on a 2 metre cross-arm. The cross-arm was attached to the top of a 7 metre tall tower. The city hall is located on a slope rising from the ocean side to an elevation of slightly higher than 80 metres above sea level over a horizontal distance of 2 kilometres. Within a radius of 1/2 kilometre, the area surrounding the site is covered with roughness elements consisting mostly of one and two storey buildings in a suburban setting. Although this site is not ideal because of the length of the tower mounted on the two storey building, it was deemed the best in the area offering the best exposure and security. The Sunset site was in the security compound of the Mainwaring substation in South Vancouver. This site is also the site used in many of the urban-meteorology/climatology studies carried out and/or supervised by Drs. Oke and Steyn of the U.B.C. Department of Geography. (Styen,1980; Steyn and Oke 1982; Cleugh and Oke, 1986) One of the reasons -185-this site was chosen for these studies was the length and quality of the fetch in all directions. The surrounding area is of very uniform roughness length and consists mainly of one and two storej' buildings, again in a suburban setting. Detailed descriptions of the site are given given by Kalanda (1979), Steyn (1980). The instruments were mounted on a 1 metre cross-arm which is on a 2 metre extension at the 30 metre height of a lattice work tower. Although the site was ideal in that it offered excellent quality and extent of fetch, there were problems with the instrumentation that caused much of the data to be spoiled. The Pitt Meadows site was located at the edge of an ultra-light aircraft facility. The site had very flat homogeneouse fetch in all directions, which consisted of either cut grass or short crops. The instrument installation has been described above. This site was of particular interest because of the location at the end of a valley occupied by Pitt Lake. There is a sharp ridge approximately 500 metres to the north of the instrument site. This ridge rises about 40 metres from the valley floor and runs parallel to the valley. There are a number of thermal and synoptic regimes that force the flow pattern within and around the Fraser Valley. The effect of these forcings should be quite evident at the end of a valley such as the one occupied by Pitt Lake. For example, one would expect to observe outflow conditions in the early mornings after a clear night on which radiative cooling has dominated. One might also expect to observe the flow of air into the interior of the province caused by a synoptic scale thermal low in the interior region. The intent of the positioning of this site was to gain more information about the effects of topographic/thermal forcing in and around-a valley such as the one containing the Pitt Meadows site. This information is of particular significance because there are a number of valleys similar to this one that are adjacent to the Fraser Valley air basin. The Aldergrove site was in a large open field adjacent to an antenna farm on the Canadian Forces Station in Aldergrove. Because of the nature of the facility, the area is large and free of obstructions that might significantly alter the flow pattern near the -186-instruments. Deciduous trees were found at the periphery of the field, more than 2 5 0 metres from the instrumentation. Few of these trees were over. 15 metres in height. The Cloverdale site was near the middle of a very large flat area (5 by 5 kilometres), south of the town of Cloverdale. The fetch was nearly ideal in all directions and consisted of short grass fields used for grazing and other agricultural activities. The tower was downed by one of the employees of the farm while mowing the field in which the instruments were located. Fortunately this happened late in the summer and did not cause prolonged interruption of data collection. The Cranberry Farm (May Brothers Farms Ltd.) site was located at the north-east end of Richmond in a bog used for growing cranberries. This particular area was quite open with the largest obstruction to the flow of air being to the north-west some 350 metres away. This obstruction consisted of a grove of trees approximately 20 metres in height and covering an area with a cross-flow length of about 60 to 70 metres. The rest of the surrounding area was exceptionally flat with very few significant obstructions to the flow. AVERAGING ALGORITHM The data from those stations using the Campbell Scientific CR21 data logger were vector averaged using the following algorithm: mean wind speed (S): N S = Z S-./N i = l 1 eq.a3.1 where Sj is the instantaneous wind speed; mean wind vector direction (0): G = t a n " 1 ( x / y ) Mean wind vector magnitude (U): U = [x2 +'y2]l/2 where: N x = Z Si sin(0j_)/N i=l N y = Z Si cos(©i)/N i=l - 1 8 8 -APPENDIX IV: RADIATION INPUT AND SLOPE CORRECTION One of the inputs needed to parameterize surface heating is the net radiation at the surface. The net all-wave radiation is arrived at by the addition of the net long and shortwave components of the radiation budget. What follows is the method used to calculate these two components. Also included in this appendix are the corrections to direct and diffuse beam radiation to account for slope angle and azimuth of sloping terrain. Longwave Radiation Component The net longwave radiation component of the radiation budget is calculated using a very simple model formulated by Idso and Jackson (1969). This model gives net radiation values from inputs of the air temperature at 10 metres above the surface. This model was chosen because only near surface temperatures are modeled. Although the model is simple, it appears to give reasonable values of net longwave radiation, (see Figure a4.1) L* = oT a 4{-0.216exp[-7.77xl0- 4(273-T a)2]} e q - a 4- 1 Shortwave Radiation Components In the interest of making the parameterization of surface heating as accurate as possible, measured values are used as a basis for solar radiation inputs. When compared to a modeled solar shortwave radiation input, this has the disadvantage of requiring more input data. The fact that the radiation terms are not modeled also makes the model less portable. The disadvantage of the increase in required input data and the decrease in portability were thought to be outweighed by the advantages of accuracy gained and the fact that, if measured data are used, the possibility exists for using the model under cloudy skies. The measured data consisted of hourly averaged values of direct normal and global solar radiation for the periods modeled. 0 --IDSO-JRCKSON ( 1 9 G 9 ) NET LONGNRVE RADIATION PARAMETERIZATION -20 -- L* vs NEAR SURFACE AIR TEMPERATURE -40 --OJ <E -60 \ 2 -80 - 1 00 20 1 40 --— t 1 1 1 1 1 1 — 2G5 270 2?5 280 285 290 295 300 TEMPERRTURE (K) Figure a4.1 Idso and Jackson (1969) Longwave Radiation Parameterization - 1 9 0 -In order to convert these data into a usable form, the values of direct beam radiation incident upon a horizontal surface are calculated. These values are then subtracted from values of global solar radiation to obtain diffuse beam solar radiation. These calculations are made using the formulation presented by Hay and McKay (1985) for the direct irradiance on an inclined surface: S = I COs(i) • eq.a4.2 c o s ( i ) = cos(ct)cos(z) + s i n ( c t ) s i n ( z ) c o s ( a - b ) eq.a4.3 cos(z) = sin(<(>)sin(6) + cos (<j>) cos (S )cos(h) eq.a4.4 cos(a) = °°s(z)sin(<fr) ~ S i n ( 5 ) eq.a4.5 sin(z)cos(<t>) Where S is the irradiance on a sloping surface, I is the radiant intensity at normal incidence, i is the angle of incidence between the sun and normal to the surface, z is the solar zenith angle, a is the slope angle,*}) is the latitude, 6 is the solar declination and h is the hour angle. As can be seen, the more complicated formulation defaults to a somewhat more simple formulation when implemented on a horizontal surface. The only inputs needed are the - 1 9 1 -solar declination (6). the latitude (a)) and the hour angle (h). The solar declination was found using a truncated Fourier series derived by Spencer (1971). Calculations of the anisotropy index ( K ) are also made at this point. As stated in equation a4.6, < is the value of the direct normal solar radiation divided by the solar constant. This anisotropy index is calculated for use in the anisotropic slope irradiance model presented by Hay and McKay ( 1985 ) . This model uses the anisotropy index in order to estimate and account for the strong anisotropy in the diffuse beam solar radiation. This model assumes that the isotropic and circumsolar components of the diffuse beam radiation are a linear combination based on the transmissivity for the direct radiation. (Hay and McKay, 1985) < = I / I n eq.a4.6 D * s = 0 . 5 D ( l + C O s ( c t ) ) ( l - < ) eq.a4.7 D S = D ( K C O S ( i ) / c o s ( z ) ) eq.a4.8 Ds = D [ > c o s ( i ) / c o s ( z ) +0.5(1.0-K)(1.0 + c o s ( o ) ] eq.a4.9 Where D g is diffuse solar irradiance for an inclined surface, D is the diffuse sky irradiance, D g is the isotropic component of the diffuse solar irradiance for an inclined surface and Dg' is is the circumsolar component of the diffuse solar irradiance for an inclined surface. - 1 9 2 -The correction represented by equations a4.2 to a4.5 is used to modify values of the direct component of the shortwave radiation, while the model shown in equations a4.6 to a4.9 is used to adjust the diffuse component for slope angle and azimuth. In order to make these data usable in the windfield model, a least squares Fourier fit is passed through the hourly data. This makes possible the interpolation of direct normal and diffuse solar radiation at any time during the modeled day. This eliminates step discontinuities in the solar heating which could occur if only the hourly values are used. Figures a4.2 through a4.11 represent measured values of direct normal, direct, global and diffuse solar radiation, as well as the anisotropy through the days of July 20 and August 23,1985. Calculated time series of direct solar radiation (horizontal surface), diffuse solar radiation and the anisotropy index are also shown. These time series are representative of a typical cloudless, midsummer day over the lower Fraser Valley. The Fourier fit through the data, as well as the coefficients used in the interpolation, are also shown. Slope Angle and Azimuth Calculation In order to make the adjustments outlined above, the values of slope angle and azimuth for all grid nodes within the domain are required. These are calculated using an algorithm devised by Steyn (1976). This method calculates the surface normal vector using a cross-product of the vectors connecting the two sets of diagonally opposed points on a grid square, (see Figure a4.13) This normal vector is then dotted with a unit vector pointing in the positive y direction to derive the zenith angle and with a horizontal unit vector pointing in the positive x direction to derive the slope azimuth. This method gives an average slope angle and azimuth for the grid square and cannot be expected to resolve variations in slope of a scale smaller than the grid square itself. JULY 20 /1985 DIRECT NORMRL SOLAR RRDIRTION OJ 800 < \ cr 600 o LO az - - - 4 - t - 7 3 2 0^ 4 00 o u 2ce F O U R I E R C O E F S . 1 = 0 , 1 0 a ( I ) x C 0 S ( 2 * P I * I / 2 4 * T ) 4 . 6 2 9 0 5 1 E + 2 - 5 . 2 2 7 0 7 5 E + 2 7 1 5 ? 1 8 E + 1 8 . 6 8 1 5 7 8 E + ! 4 . 8 0 4 9 2 1 E + 1 1 . 8 9 9 9 2 5 E + 1 7 7 1 4 9 8 E + ! 4 G 9 8 B B E - 2 2 0 5 4 B G E + 0 4 2 2 9 6 6 E + 0 •1 . 7 5 8 2 3 8 E + F O U R I E R C O E F S . 1 = 1 , 1 0 b ( I ) x S I N ( 2 * P I * I / 2 4 * T ) - 8 . 9 5 5 4 4 4 E + 1 E 2 7 8 2 1 E + 3 6 4 1 3 1 3 E + 1 7 B G 2 2 B E + 1 0 3 3 1 4 6 E + 0 4 9 8 0 1 6 E + 1 9 5 1 1 5 4 E + 0 S 3 5 7 9 4 E + 0 4 2 5 4 3 6 E + 0 4 3 4 1 9 9 E - 1 1 * MEASURED FOURIER FIT G 8 10 12 14 16 18 LOCRL RPPRRENT TIME ( H r ) < 4 t 22 21 F i g u r e a 4 . 2 D i r e c t N o r m a l S o l a r R a d i a t i o n - J u l y 2 0 , 1 9 8 5 830 - -J U L Y 2 0 / 1 9 8 5 DIRECT SOLAR R R D I R T I O N <_ 2_ \ 3 or O in u (Y 603 4 2 3 200 FOURIER COEFS. 1=0,10 a ( I ) x C0S(2*?X*I/'24*T> 2 . 839556E+2 016502E+2 0S47 15E+2 131849E+1 875 164E+0 121097E+1 020S40E-1 594151E+3 429153E-! 353914E-1 4978 SSE-1 -4 1 2 -4 -1 5 2 -2 4* * i i "i I—I I I 1 I I I I FOURIER COEFS. 1-1,10 bCI) x SIN<2xPI*I/24*T) -S.103045E+1 3.865442E+1 4 .888455E+3 -4.794596E+0 -4 . 404443E+0 -2.62S095E+0 4. 178475E+0 1 . :4S320E+0 -1 .0745G8E+3 2.689924E-1 X MEASURED FOURIER FIT S B '.3 12 14 IS 18 LOCRL RPPRRENT TIME ( H r ) 20 22 24 Fieure a4.3 Direct Solar Radiation - July 20,1985 leoa J U L Y 2 0 / 1 9 8 5 GLOBRL SOLAR R A D I A T I O N (VI < \ z o y-FOURIER COEFS. 1=0,10 a ( I ) x C 0 S ( 2 x P I * I / 2 4 * T ) 3 .31 180GE+2 .60G778E+2 . 153-43 1E+2 .556224E+1 .33 1549E+0 . 8 8 5 2 7 4 E + 0 . 3 1 0 0 7 5 E + 0 , e 4 7 5 3 4 E + 0 , 3 8 0 2 4 0 E - 1 .4G0436E -1 . 246164E -1 <J 2 0 3 - -FOURIER COEFS. 1=1,10 b ( I ) x S I N ( 2 * P I x I / 2 4 * T ) - 6 . 6 8 5 4 9 4 E + 1 3 .575027E+1 1 .0353G3E+1 - 3 . 4 8 8 2 4 1 E + 0 - 6 . 0 3 1 7 5 G E + 0 -2 .SG48G3E+0 2 .931584E+0 1.556118E+0 - 2 . 9 2 7 3 1 S E - 1 - 1 . 4 3 1 1 1 4 E - 1 U l I LOCAL APPARENT TIME (Hr) Figure a4.4 Global Solar Radiation - July 20,1985 I 323 - L 633 J U L Y 2 0 / 1 9 8 5 D I F F U S E SOLAR R A D I A T I O N a. CE _i o in UJ Ln u l - l a 630 433 233 FOURIER COEFS. 1 - 0 , 10 aCI ) x C G S ( 2 * P I * I / 2 4 * T ) 4 . 722495E+1 .372620E+1 .291061E+0 . 4786?1E+0 ,459G48E+e . 153S50E+0 . 7B05ME+0 .035187E -1 4 . 7 1 E S 3 2 E - 1 -7 . 775897E -2 - 7 . 8 G 0 8 6 1 E - 2 - 6 . 7 4 1 ! -1 - 5 . -x—*—s—X FOURIER COEFS. 1=1 ,10 b ( I ) x S I N ( 2 * P I * I / 2 4 * T ) - 6 . 5 8 i e i E E+0 - 2 .334555E+0 5 . 457B32E+0 .22269GE+0 . 6656G6E + 0 . 4 2 4 7 6 2 E - 2 172415E+0 194G 16E-1 . 5 G 3 4 5 S E - 1 .02807eE-l 1 , -1 - 7 . -1 . 4 . 7. - 4 . I X M E A S U R E D F O U R I E R F I T B 8 10 12 14 16 !B LOCRL RPPRRENT TINE (Hr) Figure a4.5 Diff u s e Solar Radiation - J u l y 20,1985 J U L Y 2 0 / 1 9 8 5 RNISOTROPY INDEX D \ rx o_ Q_ cr . 8 FOURIER COEFS. 1 = 0, 10 a ( I ) x C O S ( 2 * P I * I / 2 4 * T > 3 . 4 2 1 8 8 6 E - 1 . 884403E -1 . 6 6 7 7 6 4 E - 2 . 5 0 4 6 1 0 E - 2 . 6 1 8 2 5 1 E - 2 . 4 1 9 0 3 8 E - 2 . 3 3 1 3 4 S E - 2 .75 1437E-5 . 4 0 9 0 2 1 E - 3 . 8 1 6 0 9 1 E - 3 - 3 - 3 6 3 - t - 1 - 7 2 1 -1 .31831GE FOURIER COEFS. 1=1 ,10 b ( I ) x S I N ' C 2 * P I * I / 2 4 * T ) - 6 . 6 G 7 6 5 0 E - 2 -1 . 1 4 5 4 5 3 E - 3 2 . 7 2 8 6 ! 1 E - 2 2 . 0 3 9 5 4 8 E - 2 - 6 . 7 3 2 9 1 2 E - 3 - 1 . 8 7 G 5 0 4 E - 2 -1 . 4 9 8 3 9 2 E - 3 G . 7 0 6 2 7 2 E - 3 4 . 0 8 4 6 9 4 E - 3 1 . 8 1 6 3 6 9 E - 4 vo 6 8 1 0 1 2 14 I B I B LOCP.L APPARENT TIME (Hr) 2 0 2 2 2 4 Figure a4.6 Anisotropy Index - J u l y 20,1985 1333 CM 833 < AUG 2 3 / 1 9 8 5 DIRECT NORMAL SOLAR R A D I A T I O N S33 --\ 2 Qi cr _i o cn cr x Di 433 O Z t -u LJ Di w 233 FOURIER COEFS. 1 - 2 , 1 0 a ( I ) x C 0 S ( 2 * P I * I / 2 4 * T > 4.01S5G0E+2 - 5 . 143344E+2 2.S82178E+1 1 .004207E+2 5 .104913E+0 -3 .083049E+1 .629741E+0 .2G8188E+3 .535S41E+0 . 4 3 1 5 S 1 E - 2 . 6CS992E -1 * -t :t * FOURIER CCEFS. 1-1,10 bCI) x SINC2*?I*I/24*T) -! .0S0237E+2 2 . 55454SE+1 5.457152E+1 -8.33165SE+0 -3 .36 1959E+1 -3.6e2112E+0 ! .53S923E+1 5.581B43E+0 -5.9ie834E+0 -3 . 359795E+3 S MEASURED FOURIER F I T I t t A 'I \ I 8 23 22 24 LOCRL RPPRRENT TINE (Hr) Figure a4.7 Direct Normal Solar Radiation - August 23,1985 ieea 803 RUG 2 3/1985 DIRECT SOLRR R R D I R T I O N ru < \ 2 cr o co FOURIER COEFS. 1=0,10 a ( I ) x C 0 S ( 2 x F I * I . / 2 4 * T ) 2 .25G803E+2 - 3 . 4 4 8 3 1 8 E + 2 1 .26 1582E+2 5 .71274SE+0 -1 .6G0749E+1 - 3 . 7 3 3 3 4 0 E + 0 3 .954248E+0 1 .GGG387E + 0 - 7 . 835449E -1 - 3 . 4 8 0 8 8 8 E -- 4 . 9 7 4 1 5 4 E -FCURIER COEFS. 1=1,10 b ( I ) x S IN (2 X P I x 1 / 2 4 x T) -5 .G9722GE+1 508305E+1 2 4 4 5 7 7 E - 3 239923E+1 4 4 4 9 5 2 E - 1 753837E+0 500511E+0 G77094E+0 4450G7E+0 3 . 8 2 2 B 5 2 E - 1 vo vO I 6 8 10 12 14 16 18 LOCAL APPARENT TIME (Hr) 20 22 24 F i g u r e a 4 . 8 D i r e c t S o l a r R a d i a t i o n - A u g u s t 2 3 , 1 9 8 5 1000 --800 < \ n a: CE o <J 200 RUG 23/1985 GLOBRL SOLRR RRDIRTION FOURIER COEFS. 1 - 0 , 1 3 a ( I ) x COS(.Z*PI*l/Z**T) 2.6S3310E+2 - 3 . 9 5 2 S 2 1 E + 2 1 . 339 1 15E + 2 1.255563E+1 -1 .729127E+1 - 5 . 0 8 4 1 3 5 E + 0 2 .968512E+0 2 .433677E+0 - 7 . S 5 2 1 0 8 E - 1 - 2 . 3 0 5 6 S 5 E - 1 - 5 . 3 5 0 3 3 8 E - 1 FOURIER COEFS. 1 - 1 , 1 0 b ( I ) x S I N ( 2 * P I * I / 2 4 * T ) -5 .G193S0E+1 3.945B83E+1 5 .114840E+0 -1 .135G1BE+1 - 3 .995907E+0 4.33G4G0E+0 344G17E+0 29 1724E+0 25G75GE+3 2 . - I . -1 . 1 . 5 0 4 7 3 5 E -S 8 10 12 14 18 18 LOCAL RPPRRENT TIME (Hr) 20 ? < * < I I I I 22 24 Figure aA.9 Global Solar Radiation - August 23,1985 1003 -J-9BB < 600 --or CE _l o cn cn L . U. 203 --RUG 2 3 / 1 9 8 5 D I F F U S E SOLRR R R D I R T I O N FOURIER COEFS. 1=0, 10 a ( I ) x C 0 S ( 2 * P I * I / 2 4 K T ) 4.0G5069E+1 .518329E+1 .48! 172E+0 .921 121E+0 1 1 I 9 9 3 E - 1 - 1 . 4 0 1 9 0 9 E + 0 -9.31G11GE-1 .901072E-1 .238 1G7E-2 12749GE-1 ,440829E-2 FOURIER COEFS. 1=1,10 b ( I ) x S I N ( 2 * P I X I / 2 4 * T > - 1 . 9 8 7 2 1 9 E - 3 - 5 . 0 0 6 7 5 2 E + 0 5 .107731E+0 1.46502GE+0 - 3 . 1 B 2 9 7 3 E + 0 - 3 . 5 2 2 7 3 4 E - 1 8 . G 4 6 4 4 2 E - 1 3 . G 2 4 0 9 2 E - 1 1 . G 8 5 2 9 9 E - 1 - 2 . 2 G 5 8 0 7 E - 1 X MEASURED FOURIER FIT 6 8 10 12 14 16 18 LOCAL APPARENT TIME (Hr) Figure a4.10 Diffuse Solar Radiation - August 23,1985 RUG 2 3 / 1 9 8 5 RNISOTROPY INDEX o n \ tr CE Y FOURIER COEFS. 1-0,!0 a(I) x COS(2*PI*I/24*T) 2 37133SE -1 3 800204E -1 2 20311 IE -2 7 424919E -2 3 777037E -3 2 279532E -2 5 64-185 IE -3 5 373748E -3 2 6 14586E -3 5 49 1604E -5 5 623843E -4 FOURIER COEFS. 1-1,10 b d ) x SIN(2*PI*I/24*T) -7.838585E-2 2. 184325E-2 4.034942E-2 -6.156805E-3 -2.263S83E-2 -2.723170E-3 137845E-2 4083S8E-3 376 182E-3 484255E-3 6 S 1 3 1 2 14 I S 1 8 LOCAL APPARENT TIME (Hr) 2 0 2 2 24 Figure a 4 . 1 1 Anisotropy Index - August 23,1985 -203-Figure a4.12 Topography defined by g r i d node e l e v a t i o n s ( a f t e r Steyn,1976) - 2 0 4 -APPENDIX V: CALCULATION OF T E M P E R A T U R E TENDENCY FROM SURFACE SENSIBLE HEAT FLUX As was the case with the friction parameterization, in order to use values calculated using the parameterization scheme, they must first undergo a modification. The units of temperature tendency are units of temperature (kelvins) divided by units of time (seconds) and the units of the surface sensible heat flux are units of power (Watts) divided by units of area (metres squared). Clearly the units must match if some form of the sensible heat flux is to be used in the temperature tendency equation. To add to the complication, there must be a distribution of heat through the boundary layer which is consistent with the physics of the model and also bears some resemblance to the reality of the atmosphere. The parameterization of the temperature profile used in the model is piecewise linear from the top boundary condition to the surface. It is made up of two segments, one from the top boundary to the top of the layer of topographic influence (boundary layer) and the other runs from the top of the boundary layer to the surface. As stated in the text, the temperature structure of the upper segment is constant through time and therefore is not affected by the surface heat flux. Any variation in the vertical temperature structure and therefore the surface pressure, is restricted to the lowest layer between the top of the boundary layer and the surface. The linearity of the lower profile segment implies that any variation in temperature at the surface must be accompanied by a temperature change over the lower segment of the profile. This change over the boundary layer segment of the profile must be one that decreases linearly with height and has a value of zero at the top of the boundary layer. This also constrains the way in which the surface flux of sensible heat is spread throughout the boundary layer. Given this parameterization, the inputs of heat from the surface must be spread through the lowest layer of the atmosphere with the flux reaching zero at the upper boundary of the affected layer. Because the surface is the only level at -205-which the temperature tendency is calculated, the temperature tendenc.y at the surface must reflect the distribution of the sensible heat throughout the boundary layer. The temperature tendency within any layer can be written in terms of the vertical divergence of the sensible heat flux. This is at the exclusion of horizontal advection and any radiative flux divergence. In order to make use of this formulation, the density as a function of height must be known. This can be done using the parameterized temperature profile of the model, the hj'drostatic equation and the equation of state for a perfect gas. ^ . { p ( z ) c p T } = I eq.a5.1 T ( z ) = T S - ( Y 1 ( Z - Z S ) ) eq.a5.2 P = pRT eq.a5.3 = -pg = 3Z -pg RT(z) eq.a5.4 If. P -g dz eq.a5.5 R[T S - ( Y l ( z - Z s ) ) ] - 2 0 6 -P ( Z ) p s J s - Y l ( z - z s ) l VR v i ] eq.ao. / RT S T s Rewriting equation a5.1 using this function for density results in: eq.a5.8 At this point it is necessary to define dQ/dZ. In reality, the heat entering the atmosphere at the surface is spread evenly throughout the mixed layer by atmospheric turbulence. This vertical turbulent mixing is enhanced by the instability created by surface heating. The sensible heat flux has the effect of increasing the temperature of the mixed layer by an amount which does not vary substantially from the surface to the top of the mixed layer. This increase in temperature allows the turbulent mixing to erode the surface-based inversion and increase the height of the mixed layer. In more complex models such as those presented by Lavoie (1972) and Keyser and Anthes (1977), this parameterization of vertical temperature structure is used. These models explicitly model the height of the mixed layer and vary substantially from the one-level model in that regard. The model, used in this study has no capacity for modeling the mixed layer depth and therefore must use a much less refined parameterization for dQ/dZ. To stay consistent with the vertical temperature structure of the model, a function for dQ/dZ must be used that will act through equation a5.1 to give a temperature change that decreases linearly with elevation. Also, the temperature change must be zero at the top of the layer of topographic influence. Use of a function dQ/dZ that linearly decreases with elevation fits this requirement rather well while staying physically realistic. Using the following formulation for dQ/dZ: 9Q Q (H-z) eq.ao.9 3 Z H2 -207-and accounting for the change in density with height by the use of equation a-5.8, the temperature tendency as a function of height can be formulated as follows: 3T _ R T S T S - ^ ( z - Z s ) Q ( H - Z ) Q * ? Z ) . W t l ' -XT-As stated earlier, because the temperature at the top of the modeled boundary layer is fixed in time and the temperature below that point is parameterized by a linear temperature profile, any effects of changes in temperature throughout that layer must act through the temperature at the surface. Equation a5.10 can be evaluated at the surface to give a temperature tendenc}' at the surface that is representative of a flux spread through the boundary layer in a linear, upward decreasing fashion. 3 T S R T S Q 3t " P ^ H E Q - A 5 - U The temperature tendency as a function of elevation is plotted in Figure a5.1. The curves presented use equation 4.10 to evaluate the temperature tendency as a function of height. Upon close examination, one can see that the temperature tendency does not actually decay linearly. This is a result of the change of density of the air with elevation. Though the relation is not linear, it is very nearly linear and is used as a close approximation. 1000 -900 -800 -700 -(A o) G00 -L. 5 0 0 -(U E ~ 400 -N 300 -200 -100 -i n ELEVRTION vs MODELED d T / d t \ ^ Ts= 280.0 and 290.0 K \ y Zs= 0.0 m. >v Th= 270.0 K N< Zh= 1000 m. Nw ' T850= 265.0 K N. Z850= 1500 m. I I 1 1 1 1 1 ^ r 1 1 1 1 u I 1 1 1 1 1 1 1 1 1 1 1 3 .2 .4 .6 .8 1 ( ( d T / d t ) / Q s ) * ( 1 0 - 6 K / s e c ) Figure a5.1 Elevation vs Modeled Temperature Tendency -209-APPENDIX VI: STATISTICAL INDICES In order to validate model results, a number of statistical indices are calculated. The formulations for these indices are presented here. Root Mean Square Error (RMSE): RMSE = {(E dj2/N)}l/2 e q . a 6 . i 3 ^ where Pj is the jth value of the predicted minus the observed values, (ie. dj = Pj-Oj) Systematic Root Mean Square Error (RMSES): RMSES = {(Z I p-i - 0 l | 2 / N } l / 2 eq.a6.2 where Pjis the ordinary least squares estimate of Pj Unsystematic Root Mean Square Error (RMSEU): RMSEU (E |p\- - p,|2/ N}l/2 eq.a6.3 J = 1- J | Index of Agreement: D2 = 1 - t^ 1 ( d J - )2}/{Z i( | P J - o | + | 0 j - o | ) 2 } e q a 6 4 where o is the mean of the observed values. -210-APPENDIX VII: MODEL RUN PARAMETERS The following is a list of parameters used in both of the model runs over the real domain surface. With the exception of the time lengths, the parameters listed were also used in the exploratory modeling. Resolution east-west direction 3703.7m Resolution north-south direction ...3500.0m Universal Transmercator coordinates of the domain boundaries: north 5480000m North south 5410000m North west 570000m East east 470000m East Grid size (north-south)x(east-west) 18x27 Coefficient of diffusion of heat 2.0x10 4 m 2 s _ 1 Coefficient of diffusion of momentum 2.0x10 4m % Depth of mesoscale influence 1000m Time step 60sec Model was run on a Floating Point Systems API 64/Max Run times approx. 850sec Model code and input arrays specifying surface characteristics used in this study are available upon request from the author. - 2 1 1 -APPENDIX VIII: LIST OF SYMBOLS A- Friction multiplication factor A r Slope of standard least squares fit B r y-intercept of standard least squares fit Coefficient of drag D g Diffuse solar irradiance for an inclined surface D s* Circumsolar component of diffuse solar irradiance for an inclined surface D s Isotropic component of diffuse solar irradiance for an inclined surface D 2 Index of agreement F Friction force vector F r Froude number H Depth of layer of topographic influence I Radiant intensity at normal incidence IQ Solar constant K m Coefficient of diffusion for momentum K t Coefficient of diffusion for heat K:;; Net shortwave solar radiation L Monin-Obukhov stability parameter L:!: Net longwave radiation N Bouyancy frequency P Atmospheric pressure P r Atmospheric pressure at upper model boundary P s Atmospheric pressure at the Earth's surface R Universal gas constant S Direct solar irradiance for an inclined surface, Mean wind speed -212-Horizontal shear stress at the top of the layer of the of topographic influence Sj Instantaneous wind speed SQ Standard deviation of observed variable Sp Standard deviation of predicted variable S g Horizontal shear stress at at the Earth's surface, Irradiance on a sloping surface Q Diabatic heating Q e Latent heat flux Q n Sensible heat flux Q Net all-wave radiation T a Air temperature Air temperature at the top of the layer of topographic influence T r Air temperature at the upper model boundary T g Air temperature at the Earth's surface (10m) U Component of wind velocity in the x-direction U Mean wind vector magnitude V Component of wind velocity in the y-direction V Geostrophic wind velocity vector V J J Horizontal wind velocity vector V G Surface wind velocity vector Z n Height of the top of the layer of topographic influence Zj Depth of the boundary layer ZQ Roughness length Z r Height of the upper model boundary Zg Surface elevation c Specific heat of air -213-dj jth error f Coriolis parameter g Force due to gravity h Height of barrier, Hour angle i Angle of incidence between sun and sloping surface k von Karman constant k Unit vector in the z-direction 0j jth observation o Mean observed value 0j Standard least squares estimate of jth observed value Pj jth prediction pj Standard least squares estimate of jth predicted value r Correlation coefficient u* Friction velocity z Vertical coordinate ct Slope angle 8 Bowen ratio, Ekman turning angle Y Lapse rate within the layer of topographic influence <5 Solar declination 9 Mean wind vector direction ic Anisotropy index P Density of air 0 Transformed vertical coordinate $ Latitude i> Stability function
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- The development and application of a one-level mesoscale...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
The development and application of a one-level mesoscale atmospheric model over the Lower Fraser Valley… Ayotte, Keith Willard 1986
pdf
Page Metadata
Item Metadata
Title | The development and application of a one-level mesoscale atmospheric model over the Lower Fraser Valley of British Columbia |
Creator |
Ayotte, Keith Willard |
Publisher | University of British Columbia |
Date Issued | 1986 |
Description | This study involves the use of a one-level mesoscale numerical model to simulate thermally driven, topographically chanelled atmospheric flow over the Lower Fraser Valley of British Columbia. The model used is a modified version of a model developed by Mass Dempsey at the University of Washington. The model was modified to include friction parameterization that was spatially variable and stability dependent. The diabatic surface heating parameterization was also modified to account for spatial variability of surface thermal characteristics. Model runs of July 20,1985 and August 23,1985 are presented in this study. The model is validated using data from a measurement network over the model domain. Statistics of the validation on the two days modeled are presented. It was found that the modified model is quite sensitive to the surface thermal characteristics specified but relatively insensitive to spatial variability of surface roughness. The model is limited in its ability to resolve effects of small topographic features in high Froude number situations but the large scale flow around topographic features is well modeled. Synoptic scale thermally induced pressure gradients are found to affect ability of the model to simulate near surface flow within the domain. It is suggested that the model be nested in a lower resolution grid in order to account for this forcing. It was concluded that this modification, coupled with more accurate specification of the surface characteristics over the model domain would enhance the agreement between modeled and observed fields. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-06-19 |
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.0096676 |
URI | http://hdl.handle.net/2429/25842 |
Degree |
Master of Science - MSc |
Program |
Geography |
Affiliation |
Arts, Faculty of Geography, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1986_A6_7 A99.pdf [ 9.37MB ]
- Metadata
- JSON: 831-1.0096676.json
- JSON-LD: 831-1.0096676-ld.json
- RDF/XML (Pretty): 831-1.0096676-rdf.xml
- RDF/JSON: 831-1.0096676-rdf.json
- Turtle: 831-1.0096676-turtle.txt
- N-Triples: 831-1.0096676-rdf-ntriples.txt
- Original Record: 831-1.0096676-source.json
- Full Text
- 831-1.0096676-fulltext.txt
- Citation
- 831-1.0096676.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-0096676/manifest