The Effects of Two External Parameters onIdealized Convective Atmospheric Boundary LayerEntrainmentbyNiamh Chaparroa thesis submitted in partial fulfillmentof the requirements for the degree ofMaster of Scienceinthe faculty of graduate and postdoctoralstudies(Atmospheric Science)The University Of British Columbia(Vancouver)December 2014c© Niamh Chaparro, 2014AbstractThe atmospheric convective boundary layer has been studied for over thirtyyears in order to understand the dynamics and scaling behaviour of itsgrowth by entrainment. This enables prediction of its entrainment rate andentrainment zone depth, and so parameterizations thereof for use in globalcirculation models.Fundamentals, such as the dependence of the entrainment rate and entrain-ment zone depth on the convective Richardson number, have been estab-lished but there is still unresolved discussion about the form of these rela-tionships. Details regarding the structure of the entrainment zone continueto emerge. The variety of convective boundary layer height and entrainmentzone depth definitions adds further complexity. The study described in thisthesis aims to join this ongoing discussion.A dry, shear-free, idealized convective boundary layer in the absence of largescale winds was modeled using a large eddy simulation. The use of ten en-semble cases enabled calculation of true ensemble averages and potentialtemperature fluctuations as well as providing smooth average profiles. Arange of convective Richardson numbers was achieved by varying the twoprinciple external parameters: surface vertical heat flux and stable upperlapse rate.The gradient method for determining local convective boundary layer heightwas found to be unreliable so a multi-linear regression method was used in-stead. Distributions of the local heights thus determined were found toiinarrow with increased upper stability. Height and entrainment zone depthwere then defined based on the ensemble and horizontally averaged poten-tial temperature profile. The resulting relationships of entrainment rate andentrainment zone depth to Richardson number showed behaviour in generalagreement with theory and the results of other studies. The potential tem-perature gradient in the upper convective boundary layer and entrainmentzone was seen to depend on the upper lapse rate, as was the positive down-ward moving temperature fluctuations at the CBL top. Overall, once thesurface heat flux was accounted for by applying the CBL height as a scale,the upper lapse rate emerged as the dominant parameter influencing scaledentrainment zone depth, and potential temperature variance in the entrain-ment zone and upper convective boundary layer.iiiPrefaceThe study described in this thesis was identified and designed by myself,Niamh Chaparro. I carried out all of the model runs and analyses.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Acronyms . . . . . . . . . . . . . . . . . . . . . . . . . . . xii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Relevant Background . . . . . . . . . . . . . . . . . . . . . . . 21.2.1 The Convective Boundary Layer (CBL) . . . . . . . . 21.2.2 CBL Height . . . . . . . . . . . . . . . . . . . . . . . . 61.2.3 CBL Growth by Entrainment . . . . . . . . . . . . . . 71.2.4 The CBL Entrainment Zone (EZ) . . . . . . . . . . . 71.2.5 Modelling the CBL and EZ . . . . . . . . . . . . . . . 91.2.6 Scales and Scaling Relations of the CBL and EZ . . . 131.3 Research Questions . . . . . . . . . . . . . . . . . . . . . . . . 202 Research Approach and Tools . . . . . . . . . . . . . . . . . 222.1 Approach to Research Questions . . . . . . . . . . . . . . . . 222.2 Large Eddy Simulation (LES) . . . . . . . . . . . . . . . . . . 242.3 Handling of Output . . . . . . . . . . . . . . . . . . . . . . . 27v3 Model (LES) Evaluation . . . . . . . . . . . . . . . . . . . . . 283.1 Initialization and Spin-Up Time . . . . . . . . . . . . . . . . . 283.2 Horizontally- and Ensemble-Averaged Profiles . . . . . . . . . 323.3 FFT Energy Spectra . . . . . . . . . . . . . . . . . . . . . . . 343.4 Visualization of Structures within the Entrainment Zone . . . 353.5 Summary of Findings . . . . . . . . . . . . . . . . . . . . . . 354 Research Answers . . . . . . . . . . . . . . . . . . . . . . . . . 374.1 Entrainment Zone Structure . . . . . . . . . . . . . . . . . . . 394.1.1 Local Mixed Layer Heights (hl0) . . . . . . . . . . . . . 394.1.2 Local Turbulent Velocity and Potential TemperatureFluctuations . . . . . . . . . . . . . . . . . . . . . . . 424.1.3 Downward-Moving Warm Air at h . . . . . . . . . . . 454.1.4 Answer to Q1 . . . . . . . . . . . . . . . . . . . . . . . 484.2 Entrainment Zone Boundaries . . . . . . . . . . . . . . . . . . 494.2.1 EZ Boundaries based on the average vertical PotentialTemperature Gradient Profiles . . . . . . . . . . . . . 494.2.2 EZ Boundaries Based on Scaled Vertical Potential Tem-perature Gradient Profiles . . . . . . . . . . . . . . . . 544.2.3 EZ Boundaries Based on Scaled Heat Flux Profiles . . 564.2.4 Answer to Q2 . . . . . . . . . . . . . . . . . . . . . . . 584.3 Entrainment Rate Parameterization . . . . . . . . . . . . . . 594.3.1 Reminder of Definitions . . . . . . . . . . . . . . . . . 594.3.2 CBL Growth . . . . . . . . . . . . . . . . . . . . . . . 594.3.3 Heights Based on the Scaled Vertical Average Poten-tial Temperature Profile . . . . . . . . . . . . . . . . . 624.3.4 Heights Based on the Scaled Vertical Average HeatFlux Profile . . . . . . . . . . . . . . . . . . . . . . . . 644.3.5 Answer to Q3 . . . . . . . . . . . . . . . . . . . . . . . 645 Results in Context . . . . . . . . . . . . . . . . . . . . . . . . 665.1 Comparison of General Set-up . . . . . . . . . . . . . . . . . . 685.1.1 Significance of Grid-size . . . . . . . . . . . . . . . . . 68vi5.1.2 Horizontal Domain . . . . . . . . . . . . . . . . . . . . 695.1.3 Initial Conditions . . . . . . . . . . . . . . . . . . . . . 705.2 Entrainment Zone Structure . . . . . . . . . . . . . . . . . . . 705.2.1 Local ML Heights . . . . . . . . . . . . . . . . . . . . 705.2.2 Local ML Height Distributions . . . . . . . . . . . . . 715.2.3 Local Vertical Velocity and Potential Temperature Fluc-tuations . . . . . . . . . . . . . . . . . . . . . . . . . . 725.2.4 Downward Moving Warm Air at h . . . . . . . . . . . 725.3 Entrainment Zone Boundaries . . . . . . . . . . . . . . . . . . 745.3.1 Direct Comparison Based on the Vertical Heat FluxProfile . . . . . . . . . . . . . . . . . . . . . . . . . . . 745.3.2 General Comparison Based on the Potential Temper-ature Profile . . . . . . . . . . . . . . . . . . . . . . . 765.4 Entrainment Rate Parameterization . . . . . . . . . . . . . . 785.4.1 Direct Comparison Based on the Vertical Heat FluxProfile . . . . . . . . . . . . . . . . . . . . . . . . . . . 785.4.2 Extending Comparison to the Average Potential Tem-perature Profile . . . . . . . . . . . . . . . . . . . . . . 796 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 806.1 Major Findings . . . . . . . . . . . . . . . . . . . . . . . . . . 806.1.1 The Gradient Method for Determining Local HeightsBased on the θ Profile is Problematic. . . . . . . . . . 806.1.2 CBL Height and EZ Boundaries can be Defined Basedon the Average Potential Temperature Profile . . . . . 816.1.3 Upper Lapse-Rate is a Critical Parameter in IdealizedCBL Entrainment . . . . . . . . . . . . . . . . . . . . 816.1.4 There are Two CBL Entrainment Regimes . . . . . . 826.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . 826.2.1 Expand Local CBL Height Determination Method . . 826.2.2 Examine Resolution Effects . . . . . . . . . . . . . . . 826.2.3 Further Explore Entrainment Regimes . . . . . . . . . 836.2.4 Apply a Mass Flux Scheme . . . . . . . . . . . . . . . 83viiReferences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84A Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88A.1 Potential Temperature: θ . . . . . . . . . . . . . . . . . . . . 88A.2 Second Law of Thermodynamics . . . . . . . . . . . . . . . . 89A.3 Reynolds Decomposition and Simplification of Conservationof Enthalpy (or Entropy) for a Dry Atmosphere . . . . . . . . 89A.4 Tri-Linear Fit for Determining Local ML Height hl0 . . . . . . 90viiiList of TablesTable 2.1 Height definitions . . . . . . . . . . . . . . . . . . . . . . . 25Table 3.1 Runs in terms of w′θ′s and initial Lapse Rate γ . . . . . . . 29Table 5.1 Comparison of Grid-Sizes used in similar Studies . . . . . 68Table 5.2 Initial Conditions used in comparable LES Studies . . . . 71Table 5.3 EZ Definitions used in comparable Studies . . . . . . . . . 76ixList of FiguresFigure 1.1 Lidar Backscatter Image of the CBL . . . . . . . . . . . . 3Figure 1.2 Visualization of Entrainment from an LES . . . . . . . . . 4Figure 1.3 Idealized vertical average Profiles for a dry CBL . . . . . 5Figure 1.4 Zero-order Jump Representation of the CBL . . . . . . . 10Figure 1.5 Alternative Potential Temperature Scale for the EZ . . . 15Figure 2.1 Height Definitions . . . . . . . . . . . . . . . . . . . . . . 25Figure 3.1 Scaled Time vs Time . . . . . . . . . . . . . . . . . . . . . 30Figure 3.2 Profiles of θ, ∂θ∂z and w′θ′ . . . . . . . . . . . . . . . . . . 31Figure 3.3 w′θ′ scaled by (w′θ′)s . . . . . . . . . . . . . . . . . . . . 31Figure 3.4 θ Profiles at 2 Hours for all Runs . . . . . . . . . . . . . . 32Figure 3.5 Scaled (w′θ′)s Profiles at 3 Hours for all Runs . . . . . . . 33Figure 3.6 FFT Energy Spectra of turbulent Velocity . . . . . . . . . 34Figure 3.7 2D horizontal slices of θ′ and w′ . . . . . . . . . . . . . . 36Figure 4.1 High local ML . . . . . . . . . . . . . . . . . . . . . . . . 40Figure 4.2 Low local ML . . . . . . . . . . . . . . . . . . . . . . . . . 40Figure 4.3 Local ML Height Distributions . . . . . . . . . . . . . . . 41Figure 4.4 Two-dimensional Distributions of w′ and θ′ for all Runs . 43Figure 4.5 Two-dimensional Distributions of w′w∗ andθ′θ∗ for all Runs 44Figure 4.6 Downward-moving warm Air at h . . . . . . . . . . . . . . 45Figure 4.7 Downward turbulent velocity Fluctuations at h . . . . . . 46Figure 4.8 Positive Potential Temperature Fluctuation at h (i) . . . 46Figure 4.9 Positive Potential Temperature Fluctuation at h (ii) . . . 47xFigure 4.10 Plot of scaled EZ upper (h1h ) and lower (h0h ) Boundaries . 49Figure 4.11 ∂θ∂z Profiles with Threshold at .0002Kkm−1 . . . . . . . . 50Figure 4.12 Scaled EZ Depth (h1−h0h ) vs inverse Richardson Number(Ri−1Delta) with threshold at .0002Kkm−1 . . . . . . . . . . 51Figure 4.13 Scaled EZ depth vs inverse Richardson Number with Thresh-old at .0004Kkm−1 . . . . . . . . . . . . . . . . . . . . . 52Figure 4.14 Scaled EZ Depth vs inverse Richardson Number (Ri−1)with threshold at .0001Kkm−1 . . . . . . . . . . . . . . . 53Figure 4.15 Scaled ∂θ∂z Profiles with Threshold at .03 . . . . . . . . . . 54Figure 4.16 Scaled EZ Depth vs Ri−1 . . . . . . . . . . . . . . . . . . 55Figure 4.17 Log-log Plot of scaled EZ Depth vs Ri−1 . . . . . . . . . . 55Figure 4.18 Scaled EZ Boundaries based on the vertical Heat FluxProfile . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56Figure 4.19 Scaled EZ Depth vs Ri−1 based on the vertical Heat FluxProfile . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57Figure 4.20 Height Definitions . . . . . . . . . . . . . . . . . . . . . . 59Figure 4.21 h vs Time for all Runs . . . . . . . . . . . . . . . . . . . . 60Figure 4.22 h vs Time for all Runs on log-log Coordinates . . . . . . . 61Figure 4.23 zfh vs Time . . . . . . . . . . . . . . . . . . . . . . . . . . 61Figure 4.24 Richardson Numbers based on∂θ∂zγ . . . . . . . . . . . . . 62Figure 4.25 Scaled Entrainment Rate vs inverse Richardson Number (i) 63Figure 4.26 Richardson Numbers based on w′θ′w′θ′s. . . . . . . . . . . . 64Figure 4.27 Scaled Entrainment Rate vs inverse Richardson Number(ii) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65Figure 5.1 Illustration of EZ Potential Temperature Scale based on γ 73Figure 5.2 Plot of the Relationship between scaled EZ Depth andRichardson Number from Federovich et al.’s (2004) . . . . 75Figure 5.3 Relationship of Scaled EZ Depth to Richardson Numberfrom Brooks and Fowler’s (2012) . . . . . . . . . . . . . . 77Figure 5.4 Plots of scaled Entrainment Rate vs Richardson Numberfrom Garcia and Mellado (2014) . . . . . . . . . . . . . . 79xiList of AcronymsCBL Convective Boundary LayerDNS Direct Numerical SimulationEL Entrainment LayerEZ Entrainment ZoneFA Free AtmosphereFFT Fast Fourrier TransformGCM General Circulation ModelML Mixed LayerLES Large Eddy SimulationRi Convective Richardson Number: Ri∆ = ghθML∆θw∗2 where ∆θ = θ(h1) −θ(h0), or Riδ = ghθMLδθw∗2 where δθ = θ0(h)− θ(h)TKE Turbulence Kinetic Energyxii1. Introduction1.1 MotivationThe daytime convective atmospheric boundary layer (CBL) over land startsto grow after sunrise when the surface becomes warmer than the air aboveit. Coherent turbulent structures (thermals) begin to form and rise, sincetheir relative warmth causes them to be less dense than their surroundings,and so buoyant. The temperature profile of the residual boundary layeris neutral; i.e., potential temperature (θ, see Section A.1) increases withheight. The thermals rise to their neutral buoyancy level, overshoot, andthen overturn or recoil. Concurrently, warm stable air from the free atmo-sphere (FA) above is trapped or enveloped and subsequently mixed into thegrowing turbulent mixed layer (ML) (Stull 1988). This mixing at the top ofthe CBL is known as entrainment and the region over which it occurs, theentrainment zone (EZ, Deardorff et al. 1980). A common, simplified con-ceptual model of this case is the dry shear free CBL (Sullivan et al. 1998,Federovich et al. 2004 Brooks and Fowler 2012). This model serves as anintellectually accessible way to understand the dynamic and complex CBLand its EZ.CBL height and the prediction thereof are important for calculating theconcentration of any atmospheric species within the ML as well as the sizesof the turbulent structures. In combination with the level at which cloudscondense (lifting condensation level) knowledge of EZ depth facilitates pre-dictions pertaining to the formation of cumulus clouds. For example cloud1cover increases as more thermals rise above their lifting condensation level(Wilde et al. 1985). Parameterizations for both CBL growth and EZ depthare required in mesoscale and general circulation models (GCMs). Further-more it is an attractive goal to develop a robust set of scales for this regionanalogous to Monin-Obvukov Theory (Stull 1988, Traumner et al. 2011,Steyn et al. 1999, Nelson et al. 1989, Sorbjan 1996).Atmospheric CBL entrainment has been studied as a separate phenomenon(Nelson et al. 1989, Sullivan et al. 1998, Federovich et al. 2004, Brooks andFowler 2012) as well as within the wider topic of entrainment in geophysicalflows (Turner 1986). There is broad agreement as to the fundamental scalingparameters and relationships involved. However, the discussion as to howthe parameters are defined and measured (Brooks and Fowler 2012, Traum-ner et al. 2011) and the exact forms of the resulting relationships continues(Sullivan et al. 1998, Federovich et al. 2004 Brooks and Fowler 2012). Thisprompts me to ask the research questions I build up to in Section 1.2 andoutline in Section 1.3.1.2 Relevant Background1.2.1 The Convective Boundary Layer (CBL)The CBL grows in three stages: (i) slowly after sunrise as the nighttimeboundary layer is burned off, (ii) rapidly in the late morning as the top risesthrough the residual layer and (iii) again slowly, when the previous day’scapping inversion is reached (Nelson et al. 1989). Convective turbulence andthe dominant upward vertical motions then begin to subside as the surfacecools. While the surface is warm, buoyancy driven thermals of somewhatuniform potential temperature (θ) and tracer concentration at their coresform and entrain surrounding air laterally as they rise, as well as trappingand mixing in stable warm from above (Stull 1988, Crum et al. 1987). Un-der conditions of strong convection and weak winds, buoyantly driven tur-bulence dominates and shear-driven turbulence is insignificant (Fedorovich2and Conzemius 2001). Thermal overshoot relative to their neutral buoyancylevel, and subsequent entrainment of the warmer air from aloft augments thewarming caused by the surface vertical heat flux (w′θ′)s (see Section A.1)and results in a θ jump or inversion at the CBL top (Schmidt and Schumann1989, Turner 1986). There may also be a residual inversion from the daybefore, possibly strengthened by subsidence (Stull 1988, Sullivan et al. 1998).Lidar images such as Figure 1.1 show the overall structure of the CBL withrising thermals, impinging on the air above (Crum et al. 1987, Crum andStull 1987, Traumner et al. 2011).Figure 1.1: Lidar backscatter image of the CBL (property of Shane D.Mayor, Department of Geological and Environmental Sciences,California State University). Time is on the bottom axis, back-scatter intensity is represented by the colours red-blue. Redcorresponds to high intensity, and therefore high aerosol con-centration. Regions with negligibly low concentrations are darkblue.This has been effectively modelled using large eddy simulation (LES) bySchmidt and Schumann (1989) who used horizontal slices of turbulent po-tential temperature and vertical velocity fluctuations (θ′ , w′) at various ver-tical levels to show how the thermals form, merge and impinge at the CBL3top with concurrent peripheral downward motions. The latter is supportedin the LES visualizations of Sullivan et al. (1998). The vertical cross sectionwithin the EZ in Figure 1.2 shows the relatively cooler thermals and trappedwarmer air as well as the closely associated upward motion of cooler air anddownward motion of warmer air.Figure 1.2: Flow visualization from Sullivan et al. 1998 showing amodelled CBL thermal enveloping FA air.On average these convective turbulent structures create a fully turbulentmixed layer (ML) with eddy sizes cascading through an inertial subrangeto the molecular scales at which energy is lost via viscous dissipation (Stull1988). Here, as represented in Figure 1.3, θ is close to uniform and increaseswith respect to time due to (w′θ′)s and the downward flux of entrained sta-ble air at the inversion (w′θ′)h. ML turbulence near the ground is dominatedby warm updraughts and cool downdraughts. With proximity to the MLtop, the updraughts become relatively cool and warmer FA air from aboveis drawn downward, so in the ML w′θ′ is positive and decreasing. Directlyabove the ML the air is stable with intermittent turbulence and, on average,transitions from a uniform ML potential temperature (∂θ∂z ≈ 0) to a stablelapse rate (γ).4θML θhz∂θ0∂z = γFAEZML∆θ(a)θθ0w′θ′sw′θ′(b)Figure 1.3: Idealized vertical average profiles for a dry CBL in the ab-sence of large scale winds or subsidence. (a) θML is the averagemixed layer potential temperature. h is the height of maximumgradient in the θ profile. The initial θ profile which has a slopeγ is represented by θ0 (dotted line). The mixed layer, entrain-ment zone and free atmosphere are denoted ML, EZ and FArespectively. (b) (w′θ′)s is the surface vertical heat flux. TheEZ boundaries (dashed lines) enclose the region of negative w′θ′ .Nelson et al. (1989) outline the stages of CBL growth from when the sub-layers of the nocturnal boundary layer are entrained, until the previousday’s capping inversion is reached and a quasi-steady growth is attained.The EZ depth relative to CBL height varies throughout these stages andits relationship to scaled entrainment exhibits hysteresis. Numerical studiestypically represent this last quasi-steady phase involving a constant (w′θ′)sworking against an inversion and/or a stable γ (Schmidt and Schumann1989, Sorbjan 1996, Sullivan et al. 1998, Federovich et al. 2004, Brooks andFowler 2012, Garcia and Mellado 2014).51.2.2 CBL HeightThe ML is fully turbulent with a uniform average potential temperature (θ)which increases sharply over the EZ . Aerosol and water-vapour concen-trations decrease dramatically with transition to the stable upper FA. Soany of these characteristics can support a definition of CBL height. Nelsonet al. (1989) defined CBL height in terms of the percentage of ML air andidentified it by eye from Lidar back-scatter images. Traumner et al. (2011)compared four automated methods applied to Lidar images:• a suitable threshold value above which the air is categorized as MLair,• the point of minimum (largest negative) vertical gradient,• the point of minimum vertical gradient based on a fitted idealizedcurve,• and the maximum wavelet covariance.CBL height detection is a wide and varied field. Both Brooks and Fowler(2012) and Traumner et al. 2011 provide more thorough reviews.Numerical models produce hundreds of local horizontal points from whichsmooth averaged vertical profiles are obtained, and statistically robust re-lationships inferred. Brooks and Fowler (2012) applied a wavelet techniqueto identify the height of maximum covariance in local vertical tracer profilesin their LES study. They compared this method to the gradient method i.e.locating the height of most negative vertical tracer concentration gradient,as well as the height of minimum w′θ′ as shown later in Figure 2.1. Thislast definition is common among LES and laboratory studies where it hasbeen referred to as the inversion height (Deardorff et al. 1980, Sorbjan 1996,Federovich et al. 2004). Sullivan et al. (1998) clarified that the extrema ofthe four w′θ′ quadrants (upward warm: w′+θ′+, downward warm: w′−θ′+,upward cool: w′+θ′−, downward cool: w′−θ′−) in the EZ more or less corre-spond to the average point of maximum ∂θ∂z (see h in Figure 1.3), whereas the6point of minimum w′θ′ was consistently lower. They defined CBL height,based on local ∂θ∂z and applied horizontal averaging, as well as in two waysbased on w′θ′ for comparison.1.2.3 CBL Growth by EntrainmentThe CBL grows by trapping pockets of warm stable air between or adjacentto impinging thermal plumes. Traumner et al. (2011) summarize two cate-gories of CBL entrainment:• Non turbulent fluid can be engulfed between or in the overturning ofthermal plumes. This kind of event has been supported by the visual-izations in Sullivan et al.’s (1998) LES study as well as in Traumneret al.’s (2011) observations. In both it appeared to occur under a weakinversion or upper lapse rate (γ)• Impinging thermal plumes distort the inversion interface dragging wispsof warm stable air down at their edges or during recoil under a stronginversion or lapse rate. This type of event is supported by the findingsof both Sullivan et al. (1998) and Traumner et al. (2011).Shear induced instabilities do occur at the top of the atmospheric bound-ary layer and in some laboratory studies of turbulent boundary layers, un-der conditions of very high stability, breaking of internal waves have beenobserved. Entrainment via the former is relatively insignificant in strongconvection, and the latter has not been directly observed in real or modeledatmospheric CBLs over the range of conditions considered here (Traumneret al. 2011, Sullivan et al. 1998).1.2.4 The CBL Entrainment Zone (EZ)The ML is fully turbulent but the top is characterized by stable air withintermittent turbulence due to the higher reaching thermals. Garcia andMellado (2014) demonstrate that the EZ is subdivided in terms of length7and buoyancy scales. That is, the lower region is comprised of mostly tur-bulent air with pockets of stable warmer air that are quickly mixed, andso scales with the convective scales (see section 1.2.6). Whereas the upperregion is mostly stable apart from the impinging thermals so scaling hereis more influenced by the lapse rate (γ). In the EZ the vertical heat flux,w′θ′ , switches sign relative to that in the ML. The fast updraughts are nowrelatively cool w′+θ′−. In their analysis of the four w′θ′ quadrants Sullivanet al. (1998) concluded that the net dynamic in this region is downwardmotion of warm air (w′−θ′+) from the free atmosphere (FA) since the otherthree quadrants effectively cancel.In terms of tracer concentration, and for example based on a Lidar backscat-ter profile, there are two ways to conceptually define the EZ. It can bethought of as the range in space (or time) over which local CBL heightvaries (Crum et al. 1987) or a local region over which the concentration (orback-scatter intensity) transitions from ML to FA values (Traumner et al.2011). The latter can be estimated using either curve-fitting or wavelet tech-niques (Traumner et al. 2011, Steyn et al. 1999, Brooks and Fowler 2012).Brooks and Fowler apply a wavelet technique to tracer profiles for the de-termination of EZ boundaries, in their 2012 LES study. However, it is morecommon in numerical modelling and laboratory studies for the EZ bound-aries to be defined based on the average vertical turbulent heat flux (w′θ′)i.e. the points enclosing the negative region as shown in Figure 1.3 (Dear-dorff et al. 1980, Federovich et al. 2004, Garcia and Mellado 2014). Bulkmodels based on the representation in Figure 1.3 assume the region of neg-ative w′θ′ coincides with the region where θ transitions from the ML valueto the FA value (Deardorff 1979, Federovich et al. 2004) but no modellingstudies use the vertical θ profile to define the EZ.Since θ modeled by an LES is not strictly constant with respect to heightin the ML (Federovich et al. 2004), a threshold value for θ or its verticalgradient must be chosen to identify the lower EZ boundary. In their 20128LES study Brooks and Fowler encountered inconsistencies when determin-ing the EZ boundaries from the average tracer profile. Although their tracerprofile was quite different to a simulated CBL θ profile, this could serve ascautionary note.Our understanding of the the characteristics and dynamics of the atmo-spheric CBL entrainment zone evolves with the increasing body of measure-ment (Traumner et al. 2011, Nelson et al. 1989), laboratory (Deardorff et al.1980) and numerical studies (Deardorff 1972, Sorbjan 1996, Sullivan et al.1998, Ebert et al. 1989, Federovich et al. 2004, Brooks and Fowler 2012, Gar-cia and Mellado 2014). Parameterizations for CBL growth and EZ depthare derived based on bulk models and are evaluated using LES output andmeasurements (Federovich et al. 2004, Boers 1989). So the relationship be-tween theory, numerical simulation and measurement is inextricable and anystudy based on one must refer to at least one of the others.1.2.5 Modelling the CBL and EZBulk ModelsBulk models for the Convective Boundary layer (CBL) based on average,vertical profiles of ML quantities can be subdivided into: (i) zero-orderjump as represented in Figure 1.4 and (ii) first-(and higher) order jumpbulk models as represented in Figure 1.3. Increased order corresponds toincreasing complexity in the shape of the θ and w′θ′ profiles at the top ofthe ML.Zero-order jump bulk models assume an ML of uniform potential temper-ature (θML) topped by an infinitesimally thin layer across which there is atemperature jump (δθ) and above which is a constant lapse rate (γ). Theassumed vertical heat flux, w′θ′ , decreases linearly from the surface up,reaching a maximum negative value (w′θ′)h . This is a constant proportion9θML θhzδθ∂θ0∂z = γ(a)θ0θ−.2(w′θ′)s0 (w′θ′)sw′θ′(b)Figure 1.4: Simplified version of Figure 1.3 such that the EZ is in-finitesimely thin. (a) h is the height of the inversion and δθthe corresponding temperature jump, that is, the difference be-tween θML and θ0(h). This is different, although related, to thejump across the EZ in Figure 1.3 (∆θ). (b) The w′θ′ profile islinear and decreasing until it reaches a maximum negative valueat h of −.2(w′θ′)s. Here there is a discontinuity as it jumps tozero.of the surface value, usually −0.2(w′θ′)s (see Section 4 in Tennekes 1973 fora discussion). At the temperature inversion w′θ′ jumps to zero across theinfinitesimally thin layer. Equations for the evolution of CBL height, θMLand δθ are derived on this basis (Tennekes 1973).If the CBL height (h) is rising, air is being drawn in from the stable freeatmosphere (FA) layer above and mixed with cooler air i.e. there is an overalldecrease in enthalpy at h. The rate of decrease in enthalpy with respect totime is cpρδθ dhdt (see Section A.1) per unit of horizontal area where, in theabsence of subsidence, dhdt is the entrainment rate (we). Since the lapse rate10above the inversion is stable, Tennekes (1973) equates this enthalpy loss tothe average vertical turbulent heat flux at the inversionδθdhdt = −(w′θ′)h. (1.1)The ML warming rate is arrived at via the simplified Reynolds averagedconservation of enthalpy, for which the full derivation is shown in SectionA.3.∂θML∂t = −∂∂zw′θ′ . (1.2)Assuming w′θ′ has a constant slope this becomes∂θML∂t =(w′θ′)s − (w′θ′)hh (1.3)and the evolution of the temperature jump (δθ) depends on the rate of CBLheight (h) increase, the upper lapse rate γ and the ML warming ratedδθdt = γdhdt −dθMLdt . (1.4)An assumption about the vertical heat flux at the inversion (h), such as theentrainment ratio, closes this set(w′θ′)h(w′θ′)s= −.2 . (1.5)The relevant quantities in equations 2.2 through 2.5 are idealized, ensembleaverages. There is some variation within this class of model. For examplethe rate equation for h (entrainment relation) can alternatively be derivedbased on the turbulent kinetic energy budget (Federovich et al. 2004, Stull1976a, Stull 1976b) but they are all based on the simplified θ and w′θ′ pro-files outlined above.First (and higher) order jump models assume an EZ of finite depth at the11top of the ML, defined by two heights: the top of the ML (h0) and thepoint where FA characteristics are resumed (h1). The derivations are morecomplex and involve assumptions about the EZ i.e.:• ∆h = h1 − h0 = Constant (Betts 1974)• ∆h = h1−h0 is related to the zero-order jump at h by two right angledtriangles with opposite sides of lengths h1−h and h−h0 (Batchvarovaand Gryning 1994)• ∆h or maximum overshoot distance d ∝ w∗N where w∗ is the Deardorffconvective vertical velocity scale and N =√gθ∂θ∂z is the Brunt-Vaisalafrequency (Stull 1973)• For h0 < z < h1 θ = θML + f(z, t)∆θ where f(z, t) is a dimensionlessshape factor (Deardorff 1979, Federovich et al. 2004)Although development of these models is beyond the scope of this thesis,they are mentioned to give context to the parameterizations considered inSection 1.2.6.Numerical SimulationsNumerical simulation of the CBL is carried out by solving the Navier Stokesequations, simplified according to a suitable approximation, on a discretegrid. Types of simulations can be grouped according to the scales of motionthey resolve. In direct numerical simulations (DNS) the full range of spatialand temporal turbulence are resolved from the size of the domain down tothe smallest dissipative scales i.e. the Kolmogorov micro-scales (Kolmogorov1962). This requires a dense numerical grid and so can be computationallyprohibitive.In an LES motion on scales smaller than twice or more of the grid spac-ing are filtered out and parameterized by a sub-grid scale closure model.General circulation models (GCM) solve the Navier Stokes equations on12a spherical grid and parameterize smaller-scale processes including convec-tion and cloud cover. LES has increasingly been used to better understandthe CBL since Deardorff (1972) devised it for this purpose. Sullivan et al.(1998), Federovich et al. (2004), Ebert et al. (1989) and Brooks and Fowlerin (2012) used it to study the structure and scaling behaviour of the EZ.1.2.6 Scales and Scaling Relations of the CBL and EZLength Scale (h)Deardorff (1972) demonstrated that dominant turbulent structures in pene-trative convection scale with CBL height, which he referred to as the inver-sion height but measured as the height of minimum vertical heat flux: zfas shown later in Figure 2.1 (Deardorff et al. 1980). Since then, the distinc-tion between the two has been clarified (see Section 1.2.2) and here h refersstrictly to the height of maximum average potential temperature gradient.There are alternatives. For example turbulence based definitions, such asthe velocity variance and the distance over which velocity is correlated withitself, represent the current turbulent dynamics rather than the recent tur-bulence history as does h (Traumner et al. 2011).Deardorff (Convective) Velocity Scale (w∗)Given an average surface vertical heat flux (w′θ′)s a surface buoyancy fluxcan be defined as gθ(w′θ′)s which gives the convective velocity scale whenmultiplied by the appropriate length scale. Since the result has units m3s3 acube root is appliedw∗ =(ghθ(w′θ′)s) 13. (1.6)13Deardorff (1970) confirmed that this effectively scaled the local vertical tur-bulent velocity fluctuations (w′) in the CBL. Sorbjan’s (1996) work supportsthis, even at the CBL top. The CBL entrainment rate (we = dhdt , neglectinglarge scale subsidence) depends on the magnitude of w′ which is driven by(w′θ′)s. Stability aloft suppresses dhdt so the influence of γ is indirectly ac-counted for via h in w∗.Convective Time Scale (τ)It follows that the time a thermal, travelling at velocity scaled by w∗, takesto reach the top of the CBL i.e. travel a distance h is scaled byτ = h(ghθ(w′θ′)s) 13. (1.7)This is also referred to as the convective overturn time scale. Sullivan et al.(1998) showed a linear relationship between h and time scaled by τ . Analternative is the Brunt-Vaisala frequency i.e. the time scale associated withthe buoyant thermals overshooting and sinking (Federovich et al. 2004). Theratio of these two time-scales forms a parameter which characterizes thissystem (see Sorbjan1996 and Deardorff 1979).Temperature Scales (θ∗ and δhγ)The CBL temperature fluctuations θ′ are influenced by w′θ′ from both thesurface and the CBL top. Deardorff (1970) showed that an effective scalebased on the convective velocity scale isθ∗ = (w′θ′)sw∗ . (1.8)Whereas Sorbjan (1996) showed that with proximity to the CBL top theeffects of FA stability γ become more important. I introduce an alternative14potential temperature scale for the EZ (δhγ) in Figure 1.5. This is thedifference in the initial or background potential temperature (θ0) across theupper part of the EZ, i.e. between h and h1.θh0hh1z δhγδhθθ0Figure 1.5: Representation of the θ0 difference across the upper partof the EZ, δhγ where δh = h1−h. This serves as an alternativeto the convective potential temperature scale θ∗.Convective Richardson Number (Ri)The flux Richardson (Rf ) number expresses the balance between mechanicaland buoyant production of turbulent kinetic energy (TKE) and is obtainedfrom the ratio of these two terms in the TKE budget equation∂e∂t + U j∂e∂xj= δi3gθ(u′iθ′)− u′iu′j∂U i∂xj−∂(u′je′)∂xj− 1ρ∂(u′ip′)∂xi− (1.9)e is turbulence kinetic energy (TKE). p is pressure. ρ is density. is viscousdissipation.Rf =gθ(w′θ′)su′iu′j∂U i∂xj. (1.10)Assuming horizontal homogeneity and vertically constant subsidence yieldsRf =gθ(w′θ′)u′w′ ∂U∂z + v′w′ ∂V∂z. (1.11)15Applying first order closure to the flux terms, i.e. assuming they are pro-portional to the vertical gradients, gives the gradient Richardson number(Rg)Rg =gθ∂θ∂z(∂U∂z)2+(∂V∂z)2 , (1.12)Applying a bulk approximation to the denominator, and expressing it interms of scales yields a squared ratio of two time scalesRg =gθ∂θ∂zV ∗2L2= N2 L2V ∗2 , (1.13)where V ∗ and L∗ are appropriate velocity and length scales. In the EZbuoyancy acts to suppress buoyant production of TKE. Applying the bulkapproximation to both the numerator and denominator yields the bulkRichardson number:Rb =gθ∆θL∗V ∗2 . (1.14)A natural choice of length and velocity scales for the CBL are h and w∗giving the convective Richardson number:Ri =gθ∆θhw∗2 . (1.15)Where ∆θ can be replaced by δθ as in Federovich et al. (2004) and Garciaand Mellado (2014). Ri can also be arrived at by considering the principalforcings of the system, or from non-dimensionalizing the entrainment rela-tion derived analytically (Tennekes 1973, Deardorff 1972). It is central toany study on CBL entrainment (Sullivan et al. 1998, Federovich et al. 2004,Traumner et al. 2011, Brooks and Fowler 2012).16Relationship of Entrainment Zone Depth to Richardson NumberA relationship of the scaled entrainment zone EZ depth to Ri∆hh ∝ Rib (1.16)is arrived at by considering the deceleration of a thermal as it overshoots itsneutral buoyancy level (Nelson et al. 1989). If the velocity of the thermal isassumed to be proportional to w∗ and the decelerating force is due to thebuoyancy difference, or θ jump, then the distance the thermal overshoots(d) can be approximated byd ∝ w∗2gθML∆θ . (1.17)If the EZ depth is proportional to the overshoot distance (d) then∆hh ∝w∗2gθML∆θh = Ri−1. (1.18)Alternatively, Boers 1989 integrated the internal (U), potential (P ) andkinetic (K) energy over a hydrostatic atmosphereU = cvg∫ p00Tdp. (1.19)P = RcvU, (1.20)andK = 12∫ p00w2g dp. (1.21)p0 is the surface pressure, R and cv are the gas constant and heat capacityof dry air at constant volume. T is temperature. Initially there is a flatinfinitesimally thin inversion interface which is distorted by an impingingthermal. The resulting height difference is assumed sinusoidal and an av-erage ∆h is obtained by integrating over a wavelength. At this point, no17entrainment is assumed to have occurred and all of the initial kinetic energy(Ki) has been transferred to the change in potential energy (∆P ).Ki = Pf − Pi = ∆P (1.22)Assuming a dry adiabatic atmosphere and that the vertical velocity in thelayer below the inversion can be approximated by the Deardorff velocityscale (w∗), the following expression is reached(∆hh)2∝ T0w∗2g∆θh . (1.23)The reference temperature, T0, can be replaced by θML to give∆hh ∝ Ri− 12 (1.24)Relationship of Entrainment Rate to Richardson NumberThe relationship between scaled entrainment rate and the buoyancy Richard-son number (Ri)wew∗ ∝ Ria (1.25)is arrived at according to the zero-order jump bulk model through thermo-dynamic arguments, or by integration of the conservation of enthalpy orturbulent kinetic energy equations over the growing CBL (Tennekes 1973,Deardorff 1979, Federovich et al. 2004). It has been verified in numerouslaboratory and numerical studies (Deardorff et al. 1980, Sullivan et al. 1998,Federovich et al. 2004, Brooks and Fowler 2012), but there is still some unre-solved discussion as the the exact value of a. It seems there are two possiblevalues, −32 and −1, the first of which Turner (1986) suggested occurs at highstability when buoyant recoil of impinging thermals becomes more impor-tant than their convective overturning. Assume that an impinging thermalsupplies kinetic energy (K) per unit time and per unit area for entrainment,in terms of appropriate length and time scales L∗ and t∗ as follows18K ∝ ρL∗3U∗2L∗2t∗ , (1.26)and that the corresponding change in potential energy per unit time andarea of the rising CBL is∆P ∝ g∆ρhdhdt . (1.27)If L∗ is the penetration depth of the thermals travelling at velocity scaledby w∗ against a decelerating force g∆ρρL∗ = w∗2ρ∆ρ . (1.28)where ρ∆ρ can be replaced withθ∆θ orθδθ (see Stull (1988) page 80 and SectionA.1) and t∗ is the response time of the inversion layer to a thermal of lengthht∗ =√h θg∆θ (1.29)then assuming all of K is transferred to the change in potential energy (∆P )and using the Deardorff velocity scaled, yieldsdhdtw∗ ∝θw∗2g∆θh√θw∗2g∆θh, (1.30)i.e.wew∗ ∝ Ri− 32 . (1.31)Adding further complexity to this discussion, Federovich et al. (2004) sug-gest that a similar power law relationship (a = −1.7) can be arrived atthrough defining the θ jump across the EZ rather than at h (see Figure 1.3).191.3 Research QuestionsA simplified conceptual model of the dry, shear-free CBL in the absence oflarge scale winds is represented in Figure 1.3. The two principal externalparameters in this case, are the surface vertical heat flux (w′θ′)s and theupper lapse rate (γ) (Federovich et al. 2004, Sorbjan 1996). They haveopposing effects, that is to say (w′θ′)s drives upward turbulent velocity(w′+) and so CBL growth (we) whereas γ suppresses it. Conversely theyboth cause positive turbulent potential temperature fluctuations (θ′+) andso warming of the CBL. In the EZ the thermals from the surface are nowrelatively cool. They turn downwards as they interact with the stable FAconcurrently bringing down warmer air. Sullivan et al. (1998) demonstratedthese dynamics by partitioning w′θ′ into four quadrants. Sorbjan (1996)asserted and showed that in this region the turbulent potential temperaturefluctuations (θ′) are strongly influenced by γ whereas the turbulent verticalvelocity fluctuations (w′) are almost independent thereof. Inspired by thesetwo studies and to gain some insight into the dynamics of idealized CBLentrainment I ask Q1 (Entrainment Zone Structure): How do thedistributions of local CBL height, and the joint distributions of w′and θ′ within the EZ, vary with (w′θ′)s and γ?The relationship between scaled EZ depth and Ri∆hh ∝ Rib (1.16)has been explored and justified in field measurement, laboratory and numer-ical studies. There is disagreement with respect to its exact form, in partstemming from variation in height and θ jump definitions, but in generalits magnitude relative to h decreases with increasing Ri. Although referredto in most relevant studies and relied upon in analytical models, the ver-tical average potential temperature profile has not been used to define theEZ (Deardorff et al. 1980, Nelson et al. 1989, Federovich et al. 2004, Boers1989, Brooks and Fowler 2012). This leads me to ask Q2 (EntrainmentZone Boundaries): Can the EZ boundaries be defined based on20the θ profile and what is the relationship of the resulting depth(∆h) to Ri?A further simplification to the dry, shear-free, CBL model without large scalevelocities, is to regard the EZ depth as infinitesimely small as in Figure 1.4.The relationship of the scaled, time rate of change of h (entrainment rate:we) to Ri can be derived based on this model (Tennekes 1973, Deardorff1979, Federovich et al. 2004)wew∗ ∝ Ria. (1.25)This will be referred to as the entrainment relation. Although that there issuch a relationship is well established, discussion as to the power exponentof Ri is unresolved and results from studies justify values of both −32 and−1. See Traumner et al. (2011) for a summary and review. Turner (1986)explains this disparity in terms of entrainment mechanism such that thehigher value occurs when thermals recoil rather than overturn in responseto a stronger θ jump (or inversion). Whereas Sullivan et al. (1998) noticea deviation from the lower power (−1) at lower Ri and attribute it to theeffect of a the shape of θ within a thicker EZ. Both Federovich et al. (2004)and Garcia and Mellado (2014) show how the definition of the θ jump in-fluences the time rate of change of Ri and so effects a. Q3 (EntrainmentRate Parameterization): How does defining the θ jump based onthe vertical θ profile across the EZ as in Figure 1.3 vs at the in-version (h) as in Figure 1.4, affect the entrainment relation and inparticular a?212. Research Approach and Tools2.1 Approach to Research QuestionsGeneral SetupI modelled the dry shear free CBL and EZ using LES, specifically the cloudresolving model System for Atmospheric Modelling (SAM) to be outlinedin Section 2.2. An ensemble of 10 cases was run to obtain true ensembleaverages and turbulent potential temperature fluctuations (θ′), each casehad a domain of area 3.2 x 4.8 Km2. Grid spacing was influenced by theresolution study of Sullivan and Patton (2011) and the vertical grid withinthe EZ was of higher resolution than that applied in other comparable work.The runs were initialized with a constant (w′θ′)s acting a against a uniformγ. So, the θ jump arose from the overshoot of the thermals, rather thanbeing initially imposed as in Sullivan et al. (1998) and Brooks and Fowler(2012).Model (LES) ValuationBefore addressing the questions stated in Section 1.3 I will examine themodeled output to make sure it represents a realistic turbulent CBL inChapter 3. I will verify that the average vertical profiles are as expectedand coherent thermals are being produced. FFT energy density spectra willshow if there is adequate scale separation between the structures of greatestenergy and the grid spacing, and if realistic, isotropic turbulence is being22modelled.Entrainment Zone StructureThe EZ can be thought of in terms of the distribution of individual ther-mal heights, or local heights. Sullivan et al. (1998) measured local heightby locating the vertical point of maximum θ gradient, and observed the ef-fects of varying Ri on the resulting distributions. However this method isproblematic when gradients in the upper profile exceed that at the inversion(Brooks and Fowler 2012). Steyn et al. (1999) fitted an idealized curve to aLidar backscatter profile. This method produces a smooth curve based onthe full original profile on which a maximum can easily be located. I willapply a multi-linear regression method outlined in Vieth (2011) to the localθ profile, representing the ML, EZ and FA each with a separate line segment.From this fit, I will locate the ML top (hl0). I’ll observe how the resultingdistributions are effected by changes in (w′θ′)s and γ using histograms inSection 4.1.1Sullivan et al. (1998) broke w′θ′ into four quadrants and used this combinedwith local flow visualizations to show how CBL thermals impinge and drawdown warm air from above. Mahrt and Paumier (1984) used 2 dimensionalcontour plots of local w′ and θ′ measurements to analyze their joint dis-tributions. In his 1996 LES study Sorbjan concluded that in the EZ, θ′is strongly influenced by γ whereas w′ is practically independent thereof.Influenced by these three studies, I will use 2 dimensional histograms at hand so within the EZ to look at how the distributions of local w′ and θ′are effected by changes in γ and (w′θ′)s . I will magnify the effects of γ,by applying the convective scales, θ∗ and w∗ and focus specifically on theentrained air at h in Section 4.1.3.23Entrainment Zone BoundariesHere I define the CBL height as the location of maximum vertical θ gradientas in Figure 2.1. The lower and upper EZ boundaries are then the pointsat which ∂θ∂z significantly exceeds zero and where it resumes γ. The lowerboundary requires a choice of a threshold value which should be small, pos-itive and less than γ. Since it is somewhat arbitrary I will compare resultsbased on three different threshold values in Section 4.2. Federovich et al.(2004) and Brooks and Fowler (2012) defined the EZ in terms of the verticalw′θ′ profiles as in Figure 2.1 but disagreed on the shape of the relationshipof scaled EZ depth to Ri (equation 2.1). As well as observing this rela-tionship using the height definitions based on the θ profile, I will apply thedefinitions based on the w′θ′ profile for comparison with Brooks and Fowler(2012) and Federovich et al. (2004) in Section 4.2.3.Entrainment Rate ParameterizationAs discussed in see Section 1.2.6 the form of the entrainment relation isthought to vary based on the mechanism that initiates entrainment, whichin turn depends on the magnitude of Ri. Furthermore the ways in whichthe height and θ jump are defined have an effect. I will vary the definitionof the θ jump as outlined in Table 2.1 in order to discern between how this,and variation in initial conditions, influence the entrainment relation and inparticular a. I will reproduce this analysis using height definitions based onw′θ′ for comparison with the results of Federovich et al. (2004).2.2 Large Eddy Simulation (LES)System for Atmospheric Modelling (SAM) is a Large Eddy Simulation withcloud resolving capability (Khairoutdinov and Randall 2003). The dynami-cal framework uses the anelastic equations of motion, which in tensor nota-tion are:24θzθθ00 γ∂θ∂zh1hh0-0.2 0 1w′θ′zf0zf1zfFigure 2.1: Height definitions based on the average vertical profiles.θ0 is the initial potential temperature.Table 2.1: Definitions based on the vertical θ profile in Figure 2.1. Toobtain those based on the w′θ′ profile, replace h0, h and h0 withzf0, zf and zf1CBLHeightML θ θ Jump Rih θML = 1h∫ h0 θ(z)dz ∆θ = θ(h1)− θ(h0) Ri∆ =gθML∆θhw∗2δθ = θ0(h)− θML Riδ =gθMLδθhw∗225∂ui∂t = −1ρ∂∂xi(ρuiuj+τij)−∂∂xip′ρ +δi3B+ij3f(uj−Ugj)+(∂ui∂t)l.s.(2.1)and∂∂xiρui = 0 (2.2)The over-bar denotes the horizontal average and prime denotes fluctuationsfrom the average. B is buoyancy = −g ρ′ρ , Ug is the prescribed geostrophicwind and f is the Coriolis parameter. τij is the sub-grid scale stress tensorand the subscript l.s. denotes the prescribed large scale tendency.The prognostic thermodynamical variable is the liquid water/ice moist staticenergy (hL).∂hL∂t = −1ρ∂∂xi(ρuihL+FhLi)−1ρ∂∂z (LcPr+LsPs+LsPg)+(∂hL∂t)rad+(∂hL∂t)mic(2.3)Lc and Ls are the latent heats of condensation and sublimation. Pr, Ps andPg are precipitation fluxes of rain, snow and graupel. These terms reduceto zero in the absence of condensed water and precipitation. The subscriptsrad and mic denote tendencies due to radiation and microphysics. Theliquid/ice water static energy ishL = cpT + gz − Lc(qc + qr)− Ls(qi + qs + qg) (2.4)where qc, qr, qi, qs and qg are the mixing ratios for cloud water, rain, ice,snow and graupel. Again, these reduce to zero in the absence of condensedwater. Temperature and potential temperature are diagnosed based on thisvariable, at each time-step. A simple first-order Smagorinski closure schemeis used to parameterize the sub-grid stresses and scalar fluxes. The eddydiffusivity coefficient is based on the grid scale.26The model equations are represented discretely on a fully staggered ArakawaC-type grid which is uniform in the horizontal and stretched in the verti-cal. Integration is performed using a third-order Adams-Bashforth schemewith variable time step. Momentum is advected in flux form with secondorder differencing and conservation of kinetic energy. Prognosed scalarsare advected using a three dimensional positive definite, monotonic scheme.Lateral boundaries are periodic. The top is bounded by a rigid lid, andNewtonian damping is applied in the top third of the domain to reduce theeffects of gravity waves. Surface fluxes are computed using Monin-Obvukhovsimilarity.2.3 Handling of OutputThe model was run on parallel computers in a Linux environment usingMessage Passing Interface (MPI). 3d variable fields were output every 10 -15 minutes in binary form and converted to Network Common Data Form(NetCDF). The use of Python was enabled using the netcdf4 interface. Plot-ting was done using matplotlib (Hunter 2007). Most analyses were per-formed using NumPy and SciPy (Jones et al. 2001–). The tri-linear regres-sion method, described in Appendix A.4, for determining local ML heightwas implemented using Cython.273. Model (LES) EvaluationIn Section 3.2 vertical profiles of the ensemble and horizontally averagedpotential temperature and heat flux (θ and w′θ′) will be examined for thedevelopment of the expected three layer structure (ML, EZ and FA). Inorder to verify that there is sufficient scale separation between the energycontaining turbulent structures and the grid size and so an adequate inertialsubrange, FFT energy spectra of the turbulent velocity fluctuations will beplotted in Section 3.3. To convince the reader (and myself) that multiplecoherent thermals are being produced, 2 dimensional visualizations will beshown at the three heights: h0, h and h1 in Figure 2.1.3.1 Initialization and Spin-Up TimeAll 10 member cases of the ensemble were run on a 3.2 x 4.8 km horizon-tal domain (∆x = ∆y = 25m, nx = 128, ny = 192). Grid numbers nx,ny were chosen based on the optimal distribution across processor nodes.The vertical grid (nz = 312) was of higher resolution around the entrain-ment layer (∆z = 5m), lower below (∆z = 25m) and stretched above it(∆z = 10 to 100m). This was guided by Sullivan and Patton’s 2011 LESresolution study of the CBL that showed how grid size affects, the shapesof the average profiles in particular around the EZ, as well as the extent ofthe inertial turbulence scale sub-range. The 7 runs, summarized in Table3.1 , are all initialized with a constant surface heat flux ( (w′θ′)s ) actingagainst a uniform initial lapse rate (γ) and differ from each other based onthese two external parameters.28Table 3.1: Runs in terms of w′θ′s and initial Lapse Rate γw′θ′s / γ 10 (K/km) 5 (K/km) 2.5 (K/km)150 (W/m2) t 150/10 l 150/51100 (W/m2) t 100/10 l 100/560 (W/m2) t 60/10 l 60/5 H 60/2.5Time must be allowed for statistically steady turbulent flow to be estab-lished. Sullivan et al. (1998) recommended 10 eddy turnover times based onthe convective time scale τ = hw∗ = h(ghθML(w′θ′s)) 13, and Brooks and Fowler(2012) chose a simulated time of 2 hours. Figure 3.1 shows that for all ofthe runs, at least 10 eddy turnover times were completed by 2 simulatedhours. Each run has a distinct convective velocity scale w∗, that increaseswith time. However, dividing boundary layer height (h) by w∗ to obtain τresults in a collapse from 7 to 3 curves; one for each γ.Figure 3.2 shows that at two hours there is a measurable well mixed layer(ML) where: (i) the horizontally and ensemble averaged potential tempera-ture (θ) is constant with height, (ii) its vertical gradient ∂θ∂z is close to zeroand (iii) the vertical heat flux w′θ′ is positive and linearly decreasing withheight. Above it is an EZ where the θ profile transitions through a maxi-mum to the upper lapse rate γ and w′θ′ is negative. At three hours the EZis fully contained within the vertical region of high resolution in all runs.Figure 3.3 shows that the w′θ′ profiles are similar and are scaled well by thesurface heat flux ( (w′θ′)s ) from two hours on.1Incomplete run: EZ exceeded high resolution vertical grid after 7 hours29Figure 3.1: Plots of scaled time vs time for all runs. Scaled time isbased on the convective time scale τ and can be thought of asthe number of times an eddy has reached the top of the CBL.30Figure 3.2: Vertical profiles of the horizontally and ensemble averagedpotential temperature (θ), its vertical gradient (∂θ∂z ) and verticalheat flux (w′θ′) for the 100/5 runFigure 3.3: w′θ′ and scaled w′θ′ vs scaled height for the 100/5 run313.2 Horizontally- and Ensemble-AveragedProfilesIn Figures 3.2 and 3.4 the θ profiles exhibit an ML above which ∂θ∂z > 0 andreaches a maximum value at h before resuming γ at h1. Convective bound-ary layer CBL growth is stimulated by (w′θ′)s and inhibited by γ. In Figures3.2 and 3.5 the w′θ′ profiles decrease from the surface value, (w′θ′)s, pass-ing through zero to a minimum before increasing to zero. They are similaracross runs when scaled by (w′θ′)s. All minima are less in magnitude thanthe zero order approximation, −0.2× (w′θ′)s (Tennekes 1973), but increasewith increased γ.Figure 3.4: θ Profiles at 2 Hours for all Runs32Figure 3.5: Scaled (w′θ′)s Profiles at 3 Hours for all Runs333.3 FFT Energy SpectraIn Figure 3.6, two dimensional FFT power spectra of horizontal w′ slicestaken at three different levels (h0, h and h1 as shown in Figure 2.1) arecollapsed to one dimension by integrating around a circle of wave-numbers.Isotropy in all radial directions is assumed and k =√k2x + k2y. The resultingscalar density spectra show peaks in energy at the larger scales, cascadingthrough a substantial inertial subrange to the lower scales roughly accordingto a −53 slope, lower in the EZ. At the top of the EZ where turbulence issuppressed by stability, the slope is steeper. The peak in energy occurs atsmaller scales at h as compared to at h0, indicating a change in the size ofthe dominant turbulent structures. The spectra for the horizontal turbulentvelocity fluctuations were analogous but show lower energy as expected. Allruns produced spectra with these characteristics.Figure 3.6: Scalar FFT energy vs wavenumber (k =√k2x + k2y) forthe 60/2.5 run at 2 hours, taken at three heights (h1, h andh0). E(k) is E(kx, ky) integrated around a circle of radius k.E(kx, ky) is the total integrated energy over the 2D domain. kxand ky are number of waves per domain length.343.4 Visualization of Structures within theEntrainment ZoneHorizontal slices, at h0, h and h1 of the potential temperature and verticalvelocity fluctuations (θ′ and w′) are plotted to see the turbulent structures.Figure 3.7 shows the bottom of the EZ (h0) for the 150/10 run where co-herent areas of positive and negative temperature fluctuations correspondto areas of upward and downward moving air. In (b) and (e) the individualthermals of relatively cool air are more evident at the inversion (h) and theirlocations correspond to areas of upward motion. Most of the upward movingcool areas are adjacent to and even encircled by smaller areas of downwardmoving warm air. At h1 ((c) and (f)) peaks of cool air are associated withboth up and down-welling.3.5 Summary of FindingsEach 10 member ensemble run was allowed a period of time to develop thethree layer structure (ML, EZ and FA) as seen from the average potentialtemperature (θ) and vertical heat flux (w′θ′) profiles. The convective timescale (τ) for a thermal to reach the CBL top (h) was seen to depend on γ,signalling the importance of this external parameter. FFT spectra of turbu-lent velocity fluctuations the ML showed a satisfactory inertial subrange andseveral coherent impinging thermals were observed in the EZ at any giventime after 2 hours, indicating that realistic turbulence was being simulated.35(a)(b)(c)(d)(e)(f)Figure 3.7: θ′ (left) and w′ (right) at 2 hours at h0 (a,d), h (b,e) andh1 (d,f) for the 150/10 run. The locations of two impingingthermals are circled.4. Research AnswersSection 4.1 will focus on answering Q1 (Entrainment Zone Structure):How do the distributions of local CBL height, and the joint dis-tributions of w′ and θ′ within the EZ, vary with (w′θ′)s and γ?The distributions of local ML heights at each horizontal point, in each en-semble member, will be plotted as histograms to visualize the effects of(w′θ′)s and γ. For the same reason the joint distributions of local potentialtemperature and velocity fluctuations (θ′ and w′) at h will be plotted. Focuswill then be narrowed to the average downward moving warm quadrant at h(w′−θ′+h, w′−h where θ′ > 0 and θ′+h where w′ < 0) to examine the directinfluence of γ on entrainment.To answer Q2 (Entrainment Zone Boundaries):Can the EZ boundaries be defined based on the θ profile and whatis the relationship∆hh ∝ Rib (1.16)of the resulting scaled depth (∆hh ) to Ri?in Section 4.2, Equation 1.16 will be plotted using height definitions basedon the ∂θ∂z profile as in Figure 2.1 and Table 2.1. Since the choice of a thresh-old to determine the lower EZ boundary is somewhat arbitrary, plots will be37reproduced using two additional values. For comparison with the results ofFederovich et al. (2004) and Brooks and Fowler (2012), Equation 1.16 willbe plotted using heights based on the average heat flux (w′θ′) profile.In Section 4.3 the temperature jump will be defined in four ways to answerQ3 (Entrainment Rate Parameterization):How does defining the θ jump based on the vertical θ profile acrossthe EZ as in Figure 1.3 vs at the inversion (h) as in Figure 1.4,affect the entrainment relationwew∗ ∝ Ria (1.25)and in particular a?This analysis will involve observing how CBL height evolves in time and cul-minate in four plots representing Equation 1.25 in log-log coordinates suchthat the most suitable values of the exponent a can be identified.384.1 Entrainment Zone Structure4.1.1 Local Mixed Layer Heights (hl0)In Figures 4.1 and 4.2 the local vertical θ profiles, each at a single horizontalpoint in an individual case, exhibit a distinct ML. Above, there are sharpchanges in the profile well into the free atmosphere, due possibly to waves.These render the gradient method for determining a local CBL height, hl,unreliable. Instead a linear regression method is used, whereby three linesrepresenting: the ML, the EZ and the upper lapse rate (γ), are fit to theprofile according to the minimum residual sum of squares. Determininglocal ML height (hl0) in this way was more straight forward than the lo-cal height of maximum potential temperature gradient (hl) for the reasonsstated above.Figure 4.1 shows two local θ profiles where hl0 is relatively high. A sharp in-terface is evident indicating that this is within an active thermal impingingon the stable layer. In Figure 4.2, where hl0 is relatively low, a less definedinterface indicates a point now outside a rising thermal. Under weaker sta-bility (γ), as in Figure 4.2 (a), these inactive locations show a larger verticalregion that could be called a local EZ. In two-dimensional horizontal plots,not shown here, regions of high hl0 corresponded to regions of upward mov-ing relatively cool air at h.The distribution of hl0 represents the range over which CBL height varies inspace, so as discussed in Section 1.2.4, relates to the depth of the entrain-ment zone (EZ). Figure 4.3 (a), (b) and (c) illustrate that the distributionwidens with increasing (w′θ′)s and narrows with increasing γ. When scaledby h, the local ML height distribution narrows with increased γ, in Figure4.3 (d), (e) and (f). The upper boundary seems to be constant at about1.1(×h) , whereas the lower boundary clearly increases. When hl0s are lowerand their distribution is narrower, the scaled versions have relatively largerspacing between bins and so higher numbers in each bin. In Figure 4.3 (d),39at higher (w′θ′)s there are fewer hl0h values with higher probabilities, but thewidth of the distributions is more or less constant regardless of (w′θ′)s.(a) (b)Figure 4.1: Local vertical θ profiles with 3-line fit for the 60/2.5 (a)and 150/10 (b) runs at points where hl0 is high.(a) (b)Figure 4.2: Local vertical θ profiles with 3-line fit for the 60/2.5 (a)and 150/10 (b) runs at points where hl0 is low. The red linerepresents the ML, the blue represents the EZ and the greenrepresents the FA.40(a)(b)(c)(d)(e)(f)Figure 4.3: Histograms of local ML heights (hl0) are shown in (a), (b)and (c). Probability distributions of scaled local ML height (hl0h )are shown in (d), (e) and (f). Both sets of plots are taken at5 hours and darker shading represents higher w′θ′s. Stabilitydecreases from top to bottom i.e (a) and (d) represents runswith the highest stability (γ = 10Kkm−1).4.1.2 Local Turbulent Velocity and Potential TemperatureFluctuationsThe two-dimensional histograms of θ′ and w′ , at each horizontal point ineach ensemble case, for all runs at h are plotted in Figure 4.4 to visualizehow the distributions are influenced by changes in (w′θ′)s and γ. In orderto isolate the effects of γ, w′ and θ′ are scaled by the convective velocityand temperature scales (w∗ and θ∗) respectively and plotted in Figure 4.5.Distributions of both w′ and θ′ widen with increasing (w′θ′)s. Whereas thatof θ′ increases only slightly with increasing stability (γ) in Figure 4.4. Asexpected, γ inhibits both upward and downward w′ . The scaled version inFigure 4.5 shows damping of w′ where potential temperature fluctuationsare positive. This can be seen as the horizontal tick marks bounding the w′w∗distribution become less obscured as γ increases. Concurrently, the coolestnegative θ′θ∗ become less cool, and the warmest become warmer. So theθ′θ∗distribution shifts positively with increasing γ.Although the quadrant of overall largest magnitude is that of upward mov-ing cool air (w′+θ′−), in the EZ the net heat flux is downward moving warm(w′−θ′+) air as the other three quadrants approximately cancel. This is inline with the findings of Sullivan et al. (1998).42Figure 4.4: Two-dimensional histograms of w′ and θ′ at h for w′θ′ =150 − 60 (Wm−2) (top - bottom) and γ = 10 − 2.5(Kkm−1)(left - right) at five hours43Figure 4.5: Two-dimensional distributions of w′w∗ andθ′θ∗ at h for(w′θ′)s = 150 − 60(Wm−2) (top - bottom) and γ = 10 −2.5(Kkm−1) (left - right) at 5 hours. Tick-marks are thickenedto show the narrowing of the w′w∗ distribution whereθ′θ∗ is posi-tive, as well as the positive shift in θ′θ∗ , as γ increases.444.1.3 Downward-Moving Warm Air at hThe average downward moving warm quadrant (w′−θ′+) at h representsthe pockets of trapped or engulfed warm air that become mixed into thegrowing CBL. So its magnitude can be taken as a measure of entrainment.Figure 4.6 shows that this increases, in magnitude, in time, as well as withincreasing (w′θ′)s. Grouping according to (w′θ′)s is evident yet there is notsignificant collapse when this is applied as scale in Figure 4.6 (b). Furtherpartitioning (w′−θ′+)h into its velocity and temperature components revealsadditional complexity.2 3 4 5 6 7 8Time (h)0.140.120.100.080.060.040.020.00w′− θ′+ h (ms−1 K)100/10100/560/560/2.5150/560/10150/10(a)2 3 4 5 6 7 8Time (h)0.180.160.140.120.100.080.060.04w′− θ′+ h w′ θ′s100/10100/560/560/2.5150/560/10150/10(b)Figure 4.6: Plots of (a) the average downward moving warm air ath (w′−θ′+)h and (b) (w′−θ′+)h scaled by the average verticalturbulent heat flux (w′θ′)s vs timeFigure 4.7 shows that the velocity component w′−h where θ′h > 0, is effec-tively scaled by w∗. The curves representing θ′+h where w′h > 0 vs timedo not seem to collapse when scaled by θ∗ in Figure 4.8. Figure 4.9 showsthis component is scaled more effectively by the potential temperature scaleintroduced in Section 1.2.6, δhγ, thus indicating that the effects of γ on thepositive potential temperature fluctuations at h are more important than(w′θ′)s.452 3 4 5 6 7 8Time (h)0.070.060.050.040.030.020.01( w′−) h(whereθ′ >0) (ms−1 ) 100/10100/560/560/2.5150/560/10150/10(a)2 3 4 5 6 7 8Time (h)0.070.060.050.040.030.020.010.00( w′−) h(whereθ′ >0)w∗100/10100/560/560/2.5150/560/10150/10(b)Figure 4.7: (a) Average negative vertical turbulent velocity fluctua-tion at h (w′−h) at points where θ′ > 0 and (b) w′−h whereθ′ > 0 scaled by w∗.2 3 4 5 6 7 8Time (h)0.00.20.40.60.81.01.2θ′+ (wherew′ <0) (K) 100/10100/560/560/2.5150/560/10150/10(a)2 3 4 5 6 7 8Time (h)0.00.20.40.60.81.01.2θ′+ h(wherew′ <0)θ∗100/10100/560/560/2.5150/560/10150/10(b)Figure 4.8: (a) Average positive potential temperature fluctuation ath (θ′+h) at points where w′ < 0 and (b) θ′+h where w′ < 0scaled by θ∗.462 3 4 5 6 7 8Time (h)0.000.020.040.060.080.100.12θ′+ (wherew′ <0) (K) 100/10100/560/560/2.5150/560/10150/10(a)2 3 4 5 6 7 8Time (h)0.000.020.040.060.080.100.12θ′+ h γ(h 1−h)(wherew′ <0) 100/10100/560/560/2.5150/560/10150/10(b)Figure 4.9: (a) θ′+h at points where w′ < 0 and (b) θ′+h where w′ < 0scaled by δhγ.474.1.4 Answer to Q1Using a multi-linear regression method, the local ML heights (hl0) were de-termined. Although at each horizontal point an ML of almost uniform θbased on the local profiles is evident, the region directly above it differsdepending on location as well as from the average profile. Since there is noreliable, local definition of h, I take the distributions of local ML height (hl0)to be a measure of CBL height variance in space and so the EZ. Overall MLheight and it’s variance increase with increasing (w′θ′)s and decrease withincreasing γ. These distributions approached similarity when scaled by h,showing an increase in the lower boundary (or percentile) with increased γ.I interpret this result as an indication that increased γ results in a narrowerscaled EZ depth.Two dimensional distributions of the local turbulent fluctuations, w′ and θ′at h, widen with increasing (w′θ′)s and narrow with increased γ. Both w′and θ′ at h show some variation with γ when scaled by the convective scalesw∗ and θ∗. The distribution of w′w∗ whereθ′θ∗ is positive narrows, whileθ′θ∗shifts positively.Plots of the average downward moving warm quadrant at h ( (w′−θ′+)h ) in-dicate dependence on (w′θ′)s. Breaking (w′−θ′+)h into its two componentsreveals dependence on both (w′θ′)s and γ. The average downward movingvelocity at h ( (w′−)h ), at points where there is a positive potential temper-ature fluctuation (θ′+), show clear dependence on w∗. Whereas the averagepositive potential temperature fluctuation θ′+h where w′is scaled by, arescaled quite well by δhγ. So as one would expect, the potential temperatureof the entrained warm air depends on γ.484.2 Entrainment Zone Boundaries4.2.1 EZ Boundaries based on the average verticalPotential Temperature Gradient ProfilesIn Figure 4.10 (a) the scaled upper EZ boundaries collapse to an initialvalue of approximately 1.15, decreasing to about 1.10. The scaled lowerEZ boundaries appear grouped according to γ and increase with respect totime. So overall the scaled EZ depth (∆hh =h1h -h0h ) narrows with time.Figure 4.10: Plot of scaled EZ upper (h1h ) and lower (h0h ) boundariesbased on the average vertical potential temperature gradientprofile.The lower entrainment zone boundary h0, as illustrated in Figure 4.11 is thepoint at which the vertical ∂θ∂z profile exceeds a threshold (.0002Kkm−1),chosen such that it is positive and at least an order of magnitude smallerthan γ. As suggested by the results of Section 4.1.1, the scaled EZ depthdecreases with increasing Richardson number (Ri =gθML∆θhw∗2 as in Table2.1). However, grouping of the curves representing the relationship of scaled49EZ depth to Richardson number∆hh ∝ Rib (1.16)according to γ is evident in Figure 4.12.Figure 4.11: ∂θ∂z profiles with threshold at .0002Kkm−1. Black linesrepresent the threshold at ∂θ∂z = 0.0002 for the lower EZboundary, as well as the three lapse rates (0.0025, 0.005 and.01 Km−1 ).50Figure 4.12: Scaled EZ Depth (h1−h0h ) vs inverse Richardson Number(Ri−1Delta) with threshold at .0002Kkm−151Threshold Test for Lower EZ Boundary, h0To explore how varying the threshold value effects Equation 1.16, plots anal-ogous to Figure 4.12 were produced at two additional thresholds. In Figure4.13, a higher value (.0004Kkm−1) results in a higher h0 and so a narrowerEZ. In Figure 4.14, a lower threshold value (.0001Kkm−1) results in a lowerh0. Both of these threshold values result in grouping according to γ.Figure 4.13: Scaled EZ depth vs inverse Richardson Number withThreshold at .0004Kkm−152Figure 4.14: Scaled EZ Depth vs inverse Richardson Number (Ri−1)with threshold at .0001Kkm−1534.2.2 EZ Boundaries Based on Scaled Vertical PotentialTemperature Gradient ProfilesThe curves representing Equation 1.16∆hh ∝ Rib (1.16)collapse when the heights are defined based on the scaled vertical potentialtemperature gradient (∂θ∂zγ ) profile in Figure 4.19. Here Ri = Riδ =gθML∆θhw∗2and ∆θ = θ(h1)−θ(h0). This stems from a switch in the relative magnitudesof the vertical potential temperature gradient in the upper ML which canbe seen when Figure 4.15 is compared to Figure 4.11. So from here on allheights will be defined based on the scaled average profiles.Figure 4.15: Scaled ∂θ∂z Profiles with Threshold at .03Figure 4.17 supports an exponent b = −12 at lower values of Ri, increasingto b = −1 at higher Ri.54Figure 4.16: Plot of scaled EZ depth vs Ri−1. EZ boundaries and so∆θ = θ(h1)− θ(h0) are based on the∂θ∂zγ profile.Figure 4.17: Scaled EZ depth vs Ri−1 based on the∂θ∂zγ profile in log-log coordinates to identify likely values of the exponent b.554.2.3 EZ Boundaries Based on Scaled Heat Flux ProfilesWhen based on the vertical heat flux profile, the scaled EZ depth ( zf0−zf1zf )remains more less constant with respect to time in Figure 4.18. Figure 4.19shows little or no Ri dependence.Figure 4.18: Plot of scaled upper ( zf1zf ) and lower (zf0zf) EZ boundariesbased on vertical heat flux profile56Figure 4.19: Plots of scaled EZ depth vs Ri−1∆ . EZ boundaries and so∆θ are based on the w′θ′(w′θ′ )sprofile.574.2.4 Answer to Q2Initially, I define CBL height and EZ boundaries based on the ∂θ∂z profile. AsBrooks and Fowler (2012) point out, when using an average vertical tracerprofile there is no universal criterion for a significant gradient. So a thresholdvalue for the lower EZ boundary (h0) was chosen such that it was positive,small (i.e. an order of magnitude less than γ) and the same for all runs. Forthe sake of rigor, plots of the relationship∆hh ∝ Rib (1.16)were produced based on two additional threshold values yielding analogousresults. In all three cases curves representing Equation 1.16 grouped accord-ing to γThe importance of γ is revealed again as the curves representing equation1.16 become similar when heights are based on the scaled ∂θ∂z profile,∂θ∂zγ .Further inspection shows that this change primarily occurs at the lower EZboundary (h0) when ∂θ∂z is measured as proportion of γ. The influence of γon ∂θ∂z at h0 ties in with the influence of γ on downward moving θ′+ at hshown in Section 4.1.3. This prompts the use of the scaled profiles for theheights (h0, h, h1 and zf0, zf , zf1) in the subsequent section.These results support a varying exponent b in Equation 1.16 which is lowerin magnitude (−12) at lower Ri and approaches −1 at higher Ri. This is inline with theory and the results of comparable studies so the EZ boundarydefinitions based on the∂θ∂zγ profile are valid. For comparison with resultsfrom other studies these heights are also based on the vertical w′θ′ profilesas shown in Figure 2.1. I find no clear dependence of the scaled EZ depthon Ri within this framework.584.3 Entrainment Rate Parameterization4.3.1 Reminder of DefinitionsA key finding of Section 4.2.2 was that curves representing Equation 1.16group according to γ when heights are based on the unscaled ∂θ∂z profile andthen become similar when heights are based on∂θ∂zγ . So from here on allheights will be as in Figure 4.20 and the corresponding Richardson numbers(Ri) will be as in Table 2.1.θzθθ00 1∂θ∂z/γh1hh0-0.2 0 1w′θ′w′θ′szf0zf1zfFigure 4.20: Height definitions based on the scaled average verticalprofiles. θ0 is the initial potential temperature.4.3.2 CBL GrowthConvective boundary layer height (h) in Figure 4.21 increases, rapidly atfirst, with a steadily decreasing rate and relates to the square-root of timein Figure 4.22. Federovich et al. (2004) focus on the attainment of a quasi-steady state regime in which their zero-order model applies. Within this59regime scaled CBL height , hB−12s N32 where Bs is the surface buoyancy flux,relates to the square-root of their scaled time, tN . Over the time of theruns Bs is constant and N varies much more slowly than h. So based onFigure 4.22 I conclude that over the period during which I obtain output,all runs are in this quasi-steady state. The height of minimum vertical heatflux zf is a constant proportion of h in Figure 4.23 indicating that this pointadvances more slowly than h.Figure 4.21: h vs Time for all Runs60Figure 4.22: h vs Time for all Runs on log-log CoordinatesFigure 4.23: zfh vs Time614.3.3 Heights Based on the Scaled Vertical AveragePotential Temperature ProfileThe inverse Richardson numbers (Ri−1∆ and Ri−1δ ) in Figure 4.24 decrease intime and group according to γ. There is an overall difference in magnitudesince ∆θ > δθ.(a) (b)Figure 4.24: Inverse Richardson number vs time based on the∂θ∂zγ pro-file using ∆θ across the EZ in (a) and δθ at h in (b). See Table2.1.The entrainment rate (we = dhdt ) is determined from the slope of a secondorder polynomial fit to h(time) in Figure 4.21. When we is scaled by w∗ theresulting relationship to Ri∆wew∗ ∝ Ria∆, (4.1)plotted in log-log coordinates in Figure 4.25 (a), has an exponent a = −1 atlower Ri∆ and a = −32 at higher Ri∆.In Figure 4.25 (b) the relationshipwew∗ ∝ Riaδ (4.2)62possibly approaches a value of a = −1 at higher Riδ but a value of lowermagnitude would fit better overall.(a) (b)Figure 4.25: Scaled entrainment rate vs inverse Richardson number(Ri−1), in log-log coordinates, where Ri is based on the∂θ∂zγprofile using ∆θ across the EZ in (a) and δθ at h in (b). SeeFigure 4.20.634.3.4 Heights Based on the Scaled Vertical Average HeatFlux ProfileRichardson numbers with ∆θ and δθ based on the w′θ′ profile are compara-ble with those in Section 4.3.3 although Ri∆ shows considerable scatter inFigure 4.26 (a) than that in Figure 4.24 (a).(a) (b)Figure 4.26: Inverse Richardson number vs time based on the w′θ′w′θ′sprofile using ∆θ across the EZ in (a) and δθ at zf in (b). SeeFigure 2.1 and Table 2.1.In Figure 4.27 the axes are in log-log coordinates and all heights are basedon the scaled w′θ′ profile. The relationship of scaled entrainment rate toRi∆ in (a) shows scatter and either value of a or a value in between couldfit. Whereas the exponent in the relationship to Riδ in (b) seems to changethroughout the run(s) and a value less (in magnitude) than −1 might fitbetter.4.3.5 Answer to Q3In conclusion the relationship of scaled entrainment rate to Ri∆ based onthe∂θ∂zγ profile shows the least scatter over time and between runs in Figure4.25. Here the exponent seems to start at a value close to −1 increasingin magnitude, with higher Ri, to close to −32 . This apparent change with64(a) (b)Figure 4.27: Scaled entrainment rate vs inverse Richardson number(Ri−1), in log-log coordinates, where Ri is based on the w′θ′(w′θ′ )sprofile using ∆h across the EZ in (a) and δθ at zf in (b). SeeFigure 2.1 and Table 2.1.increased Ri mirrors that seen with Equation 1.16 in Figure 4.17. It’s pos-sible that it represents a change in entrainment mechanism as discussed inSection 1.2.6. Overall the definition of the temperature jump certainly hasan effect, ∆θ yielding a value of a higher in magnitude than δθ.655. Results in ContextMuch work has been done to develop our understanding of CBL entrain-ment, so this chapter will focus on how my results fit into the discussionestablished in the literature. I focus on six closely related publications forcomparison. Sullivan et al.’s 1998 LES study was seminal in shedding lighton CBL entrainment zone structure. Whereas Brooks and Fowler’s (2012)work contains the most recent LES results on the topic framed within an upto date review of CBL height and EZ definitions. Federovich et al. (2004)bridges LES and bulk models, while the closely related DNS study of Garciaand Mellado (2014) introduces the two-layer EZ concept and answers ques-tions regarding the scale resolution required to realistically capture CBLgrowth and EZ structure. Sullivan and Patton (2011) addressed this lastpoint using an LES, and guided the choice of grid-size in the study de-scribed in this thesis. Finally Sorbjan (1996) focused on the effects of upperlapse rate on the turbulence in the upper CBL and provided ideas uponwhich I based Section 4.1.Section 5.1 draws upon the results of Sullivan and Patton (2011) to addressthe need for high resolution in the entrainment zone (EZ) such that thesteep gradients are sufficiently represented. I present and compare those ofmy results that are pertinent and refer to how Garcia and Mellado (2014)speaks to this point. Finally, I touch upon how my domain size and initialconditions compare with those of the related LES studies and consider pos-sible implications of the similarities and differences.In Section 5.2 I describe problems encountered when using the gradient66method, as well as discuss the results obtained using my chosen method ofdetermining local ML height. All of this is set in context with the resultsof Sullivan et al. (1998) and Brooks and Fowler (2012). The influence of γon the turbulent fluctuations of vertical velocity and potential temperature(w′ and θ′) is discussed and compared with the results of Sorbjan (1996) be-fore addressing the dependence of the downward moving positive potentialtemperature fluctuations at h (θ′+ where w′ < 0) on this parameter. Anexplanation of the potential temperature fluctuation scale δhγ follows.A primary goal of this thesis was to test the average potential temperatureθ profile as a basis for defining the EZ boundaries. Before comparing theresults using heights thus defined in Section 5.3, I base all heights on thevertical heat flux (w′θ′) profile to enable direct comparison with the resultsof Brooks and Fowler (2012) and Federovich et al. (2004). I discuss similar-ities, differences and possible reasons for the latter. I then compare resultsbased on the θ profile focusing on the exponent b in Equation 1.16 and howit varies depending on Ri.Section 5.4 contains an analogous comparison to that described above. Heightsare defined, first based on the w′θ′ profile for direct comparison with theresults of Federovich et al. (2004) and Garcia and Mellado (2014) and thenbased on the θ profile. Each of these two comparisons is subdivided in orderto address the effect of, defining the θ jump across the EZ as in Figure 1.3,vs at h as in Figure 1.4. In all there are four plots of Equation 1.25 toshow how the exponent a varies depending, on θ jump definition, as wellas Ri. This variation is discussed in the context of the results of the othercomparable studies.675.1 Comparison of General Set-up5.1.1 Significance of Grid-sizeSullivan and Patton (2011) found that the shapes of the average potentialtemperature (θ) and heat flux (w′θ′) profiles, as well as the measured CBLheight vary depending on grid size. The resolution at which convergencebegins is listed in Table 5.1. At lower resolution the θ and w′θ′ profiles aresuch that the entrainment zone (EZ) is a larger portion of the CBL andmeasured CBL height is higher. Overall they concluded that vertical res-olution was critical. This compliments the conclusion Brooks and Fowler(2012) reached when discussing their resolution test. That is, to capture thesteep vertical gradients in the EZ requires high vertical resolution.Table 5.1: Grid spacing around the EL used in comparable LES stud-ies. Those used for resolution tests are not listed here. For Sul-livan and Patton’s 2011 resolution study I list the grid sizes atwhich profiles within the EL and CBL height evolution began toconverge.Publication ∆x, ∆y, ∆z Horizontalin the EZ (m) Domain (km2)Sullivan et al. (1998) 33, 33, 10 5 x 5Federovich et al. (2004) 100, 100, 20 5 x 5Brooks and Fowler (2012) 50, 50, 12 5 x 5Sullivan and Patton (2011) 20, 20, 8 5 x 5This study 25, 25, 5 3.4 x 4.8As Turner discusses in his 1986 review of turbulent entrainment, smallerscale processes such as those at the molecular level are relatively unimpor-tant. Large scale engulfment and trapping between thermals dominates. Ifthe ergodic assumption holds and potential temperature variance (θ′2) iscalculated based on the difference at a point from the horizontal average,68it is a measure of horizontal variance. Sullivan and Patton (2011) foundthat the vertical distance over which θ′2 varied significantly, more or lessconverged at the resolution shown in Table 5.1. However, the maximumθ′2 continued to increase up to their finest grid spacing (∆x = 5, ∆y = 5,∆z = 2).The question as to whether mixing and gradients within the EZ are ade-quately resolved motivates DNS studies such as that of Garcia and Mellado(2014). These authors found the entrainment ratio (w′θ′zfw′θ′s) to be about 0.1which is lower than that observed by Federovich et al. (2004), but close towhat was seen here in Figures 3.5 and 3.2. Based on their w′θ′ profilesthe depth of the region of negative flux is comparable to what is shownin Figure 4.18. Furthermore, these authors concluded that the productionand destruction rates of turbulence kinetic energy (TKE), as well as theentrainment ratio used to calculate the entrainment rate, were effectivelyindependent of molecular scale processes.The FFT energy spectra (Figure 3.6) of the turbulent velocities at the topof the ML show a substantial resolved inertial subrange giving confidence inthe choice of horizontal grid size used. In the EZ where turbulence is inter-mittent, the dominant energy containing structures are smaller, and decayto the smallest resolved turbulent structures is steeper. This is consistentwith the assertion of Garcia and Mellado (2014) that the EZ is separatedinto two sub-layers in terms of turbulence scales.5.1.2 Horizontal DomainThe horizontal domain in this study is relatively small (see Table 5.1). How-ever, visualizations of horizontal and vertical slices clearly showed multi-ple resolved thermals. Their diameters increased with CBL height, butremained less than or on the order of the height of the CBL. Sullivan et al.(1998) carried out one run on a smaller domain with higher resolution, no-69ticed it resulted in lower CBL height and concluded this was due to restrictedhorizontal thermal size. However, given the results of Sullivan and Patton(2011) discussed in Section 5.1.1 it could have been an effect of grid-size.When defining heights based on average profiles Sullivan et al. (1998) pro-duced jagged, oscillating time-series and Brooks and Fowler (2012) encoun-tered significant scatter in plots of Equation 1.25. But the heights based onaverage profiles here, using an ensemble of cases, varied smoothly in timein Figure 4.21. This could be attributed to a smoother profile based on agreater number horizontal points (10*128*192).5.1.3 Initial ConditionsThe principle parameter describing the balance of forces in dry, idealizedCBL entrainment is the Richardson number (Ri) and its magnitude de-pends on the way in which the θ jump is defined. Varying the θ jumpdefinition causes identical conditions to be described by different Ri values.The Ri range in this study was dependent on variation in γ (see Figure 4.24).Brooks and Fowler (2012) and Sullivan et al. (1998) imposed a θ jump ofvarying strength topped by a constant γ. Whereas Federovich et al. (2004)initialized with a layer of uniform θ. They varied γ and kept w′θ′s constantfor each run. Their initial conditions, definitions of the θ jump and Ri rangeare directly comparable to those of this study, whereas those of Brooks andFowler (2012) and Sullivan et al. (1998) are quite different.5.2 Entrainment Zone Structure5.2.1 Local ML HeightsSullivan et al. (1998) determined local CBL height by locating the point ofmaximum ∂θ∂z . Analysis of the resulting distributions showed dependence ofstandard deviation and skewness on Richardson number (Ri). The normal-70Table 5.2: Initial Conditions used in comparable LES StudiesPublication w′θ′s γ Initial θ RiWm−2 Kkm−1 Jump K rangeSullivan et al. (1998) 20 - 450 3 .436 - 5.17 1 - 100Federovich et al. (2004) 300 1 - 10 NA 10 - 40Brooks and Fowler(2012)10 -100 3 1 - 10 10 - 100This study 60 - 150 2.5 - 10 NA 10 - 30ized standard deviation decreased with increased Ri whereas skewness wasalmost bimodal; being negative at high Ri and positive and low Ri. Initiallyin this study, I applied a similar method and found local CBL height distri-butions with lower Ri to have positive skew. Upon exhaustive inspection oflocal vertical θ profiles such as those in Figure 4.2, it became evident thatat certain horizontal points high gradients well into the free atmosphere ex-ceeded those closer to the location of the CBL height reasonably identifiedby eye.5.2.2 Local ML Height DistributionsLocating the local ML height (hl0) using the multi-linear regression methoddescribed in Chapter 2 proved more reliable than the gradient method dis-cussed above. The resulting distributions, normalized by CBL height, h,(see Figure 2.1) in Figure 4.3, showed a decrease in the lowest hl0h resultingin an apparent increased negative skew with decreasing stability (decreas-ing Ri). This, combined with a widening of the distribution agrees, withthe findings of Sullivan et al. (1998) and supports the results based on theaverage profiles in Section 4.2. The approximate scaled EZ based on the hl0hdistributions is about 0.2 - 0.4 whereas that based on distributions of localmaximum tracer gradients by Brooks and Fowler (2012) was smaller (0.05- 0.2). However, the local maximum gradient of the tracer profile would71likely be within the EZ at points outside an actively impinging plume andso higher than hl0 defined here.5.2.3 Local Vertical Velocity and Potential TemperatureFluctuationsAs expected, with increased w′θ′s the variance and magnitude of the verticalvelocity fluctuations within and at the limits of the EZ increase. Greaterturbulent velocity causes a higher CBL and a deeper EZ over which, rela-tively warmer air from higher up is brought down, and relatively cooler airfrom below is brought up. So the magnitude of the potential temperaturefluctuations (θ′) and the width of their distribution increases. All of thisagrees with the findings of Sorbjan (1996), but the portion of the scaled w′(w′w∗ ) distribution where scaled θ′( θ′θ∗ ) is positive, in Figure 4.5, appears tonarrow as γ increases. This is somewhat at odds with his assertion that w′is independent of this parameter while the effectiveness of w∗ as a scale forw′− where θ′ > 0 in Figure 4.6 supports it.5.2.4 Downward Moving Warm Air at hAlthough the motion of the thermals dominates within the EZ, the w′−θ′−,w′+θ′− and w′−θ′− quadrants do approximately cancel leaving w′−θ′+ asthe net dynamic, as Sullivan et al. (1998) concluded. The downward mov-ing warm quadrant at h ( (w′−θ′+)h ) represents warmer free atmosphere(FA) air that is being entrained. So its magnitude, at a certain point intime, is an indication of how much the region below will be warmed dueto entrainment at a successive time. The increase of (w′−θ′+)h in time isprimarily due to the increased average positive potential temperature fluc-tuation at h ( (θ′+)h ) which is effectively scaled by the temperature scale(h1 − h)γ = δhγ (see the Figures of Section 4.1.3). A similar scale wasintroduced by Garcia and Mellado (2014) to further their line of reasoningthat the buoyancy in the upper EZ is determined by γ. Figure 5.1 illustrates72a broad qualitative explanation. At h much of the air is at the background(or initial) potential temperature θ0(h). Some air at potential temperatureθ = θ0(h1) is brought down from the upper EZ limit (h1) resulting in posi-tive potential temperature fluctuations (θ′+) at h.Figure 5.1: Illustration of the potential temperature scale (h1−h)γ =δhγ: The curves represent a vertical cross-section of thermaltops. Between them is stable air at the initial lapse rate γ.h1 and h correspond approximately to the highest and aver-age thermal height respectively and h0 is the top of the wellmixed region (ML). The horizontally uniform, initial potentialtemperature is θ0 = θ0. A thermal will initiate the downwardmovement of air from h1 to h, and the difference between itspotential temperature and that of the background stable air ath is (h1 − h)γ = δhγ.Garcia and Mellado (2014) suggest that the buoyancy in the lower portionof the EZ, i.e. from a point just below h down, is more strongly influencedby the vigorous turbulence of the ML than by γ. So mixing reduces thedifference between, the potential temperature at the top of the ML, andthat at or just below h. However the observation in Section 4.2.2, that themagnitude of the average vertical potential temperature gradient (∂θ∂z ) in theupper ML increases with increasing γ, indicates that the influence of thisparameter extends further down. On a related note, the magnitude of theminimum heat flux ( (w′θ′)zf ) is seen to increase with increasing γ, here73and in both Sorbjan (1996) and Federovich et al. (2004). It is reasonableto suggest this leads to an increased negative vertical heat flux gradient(−∂w′θ′∂z ) in the lower EZ and so increased warming per Equation A.21.∂θ∂t = −∂∂zw′θ′ (A.21)5.3 Entrainment Zone BoundariesThe EZ is inhomogeneous, but on average is a region of transition as clearlyrepresented by the θ profile. It’s where relatively cooler thermals overturn orrecoil initiating entrainment as represented by the vertical heat flux (w′θ′)profile. The θ profile partially characterizes the thermodynamic state of theCBL as well defining its three layer structure. It is directly comparable toboth bulk models and local θ profiles which in turn are comparable to asounding, unlike a w′θ′ profile which is an inherently average quantity.5.3.1 Direct Comparison Based on the Vertical Heat FluxProfileNeither of the two comparable LES studies in Table 3.3 define the EL basedon the ∂θ∂z profile. So to enable direct comparison, heights were based onthe heat flux (w′θ′) profile as in Figure 4.20. In this framework Federovichet al. (2004) show decreasing scaled EZ with increasing Ri and conclude anexponent b = −12 . They attribute the decrease in the overall scaled depth toa slight decrease in the scaled upper boundary over time. However based ontheir plot in Figure 5.2 the decrease seems more than slight, varying fromabout 0.5 to 0.2.Brooks and Fowler (2012) found no clear Ri dependence of the scaled EZdepth defined based on the w′θ′ profile. But their definition hinged solelyupon the lower part (zf1 − zf ) which according to Federovich et al. (2004)does not vary in time. Figure 4.19 of this thesis shows that when I definedthe EZ based on the w′θ′ profile as Federovich et al. (2004) did, the scaled74Figure 5.2: Figure 9 from Federovich et al. (2004) representing Equa-tion 1.16 using three different Richardson numbers, in log-logcoordinates. Heights are based on the w′θ′ profile as in Figure4.20 and their zi is my zf . δzizi is then the scaled EZ depth.Ri∆b (circles) and Riδb (crosses) correspond directly to thosedetermined here using δθ and ∆θ. Note that their ∆ refers tothe smaller jump measured at zf , whereas I use it for the larger.RiN (triangles) is the Richardson number defined in Equation1.13, with w∗ and zf as the velocity and length scale.EZ depth had no clear dependence on Ri. This is supported by the similar-ity in time and across runs of the vertical heat flux profiles when scaled by(w′θ′)s in Figures 3.3 and 3.5.The most obvious possible cause for disagreement with the results of Federovichet al. (2004) is the difference in grid size shown in Table 5.1. Inspection oftheir w′θ′ profiles confirms a relatively deeper scaled region of negative fluxas compared with those seen here ( .4 vs .25). Their surface heat flux w′θ′swas twice the highest used here, but their range of Ri is comparable to thatof this study. The latter point although not directly relevant here, serves asconfirmation that γ is the more influential parameter.75Table 5.3: EZ Definitions used in comparable Studies. See Figure 4.20Publication EZ Depth CBLheightθ JumpFederovich et al. (2004) zf1 − zf0 zf θ(zf1)− θ(zf0)Brooks and Fowler(2012)2×(zf−zf0) zf average of localvalues5.3.2 General Comparison Based on the PotentialTemperature ProfileHere, when heights are defined based on the scaled vertical potential tem-perature gradient profile∂θ∂zγ the curve representing Equation 1.16∆hh ∝ Rib (1.16)shows an exponent b which increases in magnitude, from about −12 as pre-dicted and seen by Boers (1989), to about −1 as justified in Nelson et al.(1989), with increasing Ri (decreasing Ri−1). Overall there is a clear nar-rowing of the scaled EZ depth with increasing Ri (decreasing Ri−1) as sup-ported by the local height distributions in Section 4.1.1. Although based ondifferent height definitions, Federovich et al. (2004) concluded an exponentb = −12 and Brooks and Fowler’s (2012) plots show curves with an apparentexponent less in magnitude than −1, in Figure 5.3.The curves representing each run in Figure 5.3 fan out. In Figure 4.12 of thisthesis, before scaling, the ∂θ∂z profile curves separate out, but in the reverseorder. CBLs under higher stability, and so higher Ri, have larger scaled EZdepths. Whereas Brooks and Fowler’s (2012) runs with initially lower Rihave larger scaled EZ depths than those with higher, even where Ri valuesoverlap. Nonetheless, that there appears a family of separate but similar76Figure 5.3: Panel (a) from Figure 5 in Brooks and Fowler (2012) rep-resenting Equation 1.16:The normalized EZ depth is determinedin three ways (i) the upper and lower percentiles from the dis-tribution of local CBL height (maximum tracer gradient), nor-malized by the average of the local heights (pale grey) (ii) theaverage of local scaled EZ depths based on wavelet covariance(dark grey) and (iii) the average of the locally determined EZdepths scaled by the average of the locally determined heights(black), based on wavelet covariance. Their θ jump is an aver-age of the potential temperature differences across the local EZdepths.curves rather than a single curve hints at an underlying scaling parameter.Neither study referenced in Table 3.3 addresses the change in exponent withincreased Ri that I observe in Figure 4.17. It is reasonable to suggest thatthis represents a change in entrainment mechanism. Sullivan et al. (1998)observed enfolding and engulfment at lower Ri. Whereas at higher Ri whenmotion is more restricted, entrainment seemed to occur via trapping of thin-ner wisps at the edge of an upward moving thermal. Turner (1986) also dis-tinguishes between entrainment by convective overturning and recoil. Garciaand Mellado (2014) refer to a change in entrainment rate due to the effects77of increased stability on the upper EZ sub-layer. In this study, the nar-rowing of the EZ depends predominantly on the magnitude of the averagevertical potential temperature gradient ∂θ∂z in the lower EZ and upper ML.However, the scaled magnitude of upper limit in Figure 4.10 does appearto decrease slightly in time. This could correspond to the slowly decreasingupper sub layer of the EZ mentioned in both Garcia and Mellado (2014)and Federovich et al. (2004).5.4 Entrainment Rate ParameterizationRi magnitude determined in this and the comparable studies is primarilyinfluenced by the magnitude of the θ jump. Here, I define it in two ways asFederovich et al. (2004) did. I do this based on the w′θ′ profile, as in Figure4.20 and Table 2.1 for the purpose of direct comparison and to observe howthe change in definition effects Equation 1.25.wew∗ ∝ Ria (1.25)5.4.1 Direct Comparison Based on the Vertical Heat FluxProfileThe larger jump, i.e. that taken across the EZ (∆θ) as in Figure 1.3, yieldsa larger value of a as Federovich et al. (2004) conclude. Garcia and Mel-lado (2014) interpret both curves as asymptotic to straight lines (a = −1)as the upper EZ sub-layer narrows. Based on their plots in Figure 5.4, inthe absence of their justification based on the derivation of the entrainmentrelation, for ∆θ I see a curve (grey and blue) with increasing exponent ex-ceeding magnitude −1 at higher Ri. For their version of δθ I see a curve(grey and red) with exponent less in magnitude than −1.78Figure 5.4: Figure 11 from Garcia and Mellado (2014) and represent-ing equation 1.25 based on the two θ jumps. The grey and bluecurve is based on ∆θ and the (grey and) red curve is based onθ(h)−θ0(h) which is slightly different to the δθ defined here andin Federovich et al. (2004). The dashed and continuous blacklines represent the straight lines to which the curves asymptoteaccording to their analysis. Their heights are comparable tothose based on the heat flux (w′θ′) profile in Figure 4.20.5.4.2 Extending Comparison to the Average PotentialTemperature ProfileThere is an analogous distinction between curves representing Equation 1.25using ∆θ and those using δθ, when all heights are based on the∂θ∂zγ profile.Scatter is least when the θ jump is defined across the EZ. In Figure 4.25a = −32 fits at higher Ri (lower Ri−1) and a = −1 seems to fit at lower Ri.Combined with the apparent change in b for Equation 1.16 I interpret thisas an indication of a change in entrainment regime at increasing Ri.796. Conclusion6.1 Major FindingsThe dry idealized convective atmospheric boundary layer (CBL) was mod-eled using large eddy simulation (LES). Although this has been done beforeand a broad understanding of the dynamics and scaling behaviour has beenestablished, discussion of details continues. This study was intended to con-tribute to this discussion and shed light on some of these details. It wasguided by the questions outlined in Section 1.3 and answered in Chapter 4and concludes with the following four points:6.1.1 The Gradient Method for Determining Local HeightsBased on the θ Profile is Problematic.Local θ profiles vary depending on location. The top of an active thermalimpinging on the free atmosphere (FA) as in Figure 4.1 is characterizedby a steep gradient comparable to the zero-order model representation inFigure 1.4. At other locations, for example where a thermal has overturnedor recoiled and some entrainment has been initiated as in Figure 4.2, thereis a region over which the θ profile transitions to the upper lapse rate (γ).That is, there is a local entrainment zone (EZ). At such locations, thereare gradients well into the FA that exceed any within the EZ, as well as anabsence of a well-defined local CBL height. This presents both a practicaland conceptual challenge to the gradient method, while determination ofthe ML using piecewise linear regression is more reliable.806.1.2 CBL Height and EZ Boundaries can be DefinedBased on the Average Potential Temperature ProfileThe θ profile characterizes the dry, idealized CBL and links bulk models tosoundings via an LES. Both the EZ depth and CBL height based on theaverage∂θ∂zγ profile showed dependence on Ri (Sections 4.2 and 4.3) as seenin other studies and justified theoretically. So this is a valid way of definingthe CBL and its EZ.6.1.3 Upper Lapse-Rate is a Critical Parameter inIdealized CBL EntrainmentThe magnitude and variance of local ML height, increase with increasingw′θ′s, and decrease with increasing γ. The same can be said for the ver-tical velocity fluctuations (w′) in the EZ. However, increased γ results inan increase in the positive potential temperature fluctuations (θ′+) at h.The magnitude of (θ′+) at points where w′ is negative represents downwardmoving entrained air and depends on γ (Section 4.1). Below h, in the lowerEZ, the average vertical potential temperature gradient (∂θ∂z ) also dependson γ (Section 4.2.2). So, the growth of the idealized dry CBL is driven by(w′θ′)s and suppressed by stability (γ). But CBL warming is due, in part,to the entrainment of air from aloft the potential temperature of which inturn depends on γ.Evidence for the influence of γ is seen throughout this study. Distribu-tions of scaled local ML heights approach similarity, when γ is constant but(w′θ′)s is varied (Figure 4.3). Curves representing Equation 1.16 group ac-cording to γ when based on the ∂θ∂z profile, but collapse once based on∂θ∂zγ(Section 4.2.2). The convective time scale τ = w∗h and Ri group accordingto γ (Figure 3.1) lending support to Federovich et al.’s (2004) use of theBrunt-Vaisala time scale. It seems that once the effect of the surface heatflux (w′θ′s) is accounted for through h, γ emerges as the dominant param-eter in dry, idealized CBL entrainment.816.1.4 There are Two CBL Entrainment RegimesTurner (1986) outlined and theoretically justified two distinct convectiveboundary layer entrainment regimes wherein the scaled entrainment rateshave different Ri dependence. The LES flow visualizations of Sullivan et al.(1998) showed large scale engulfment at lower Ri. At higher Ri, trapping ofsmaller volumes of stable air between and at the edges of impinging thermalsappeared to be the dominant mechanism. The CBL entrainment zone mea-surements analyzed in Traumner et al. (2011) further support the conceptof varying entrainment mechanism depending on the strength of the upperlapse rate γ. Finally, both Federovich et al. (2004) and Garcia and Mellado(2014) discuss the varying dependence of the scaled entrainment rate on Rias the effects of upper stability become more important. On these groundsI attribute the change in exponent in the plots of equations 1.16 and 1.25 inFigures 4.17 and 4.25 to a change in entrainment regime as Ri increases.6.2 Future WorkSome ideas as to how the work in this thesis could be completed or extendedare as follows:6.2.1 Expand Local CBL Height Determination MethodTo further the tri-linear regression method described in Section A.4, the EZcould be approximated by a suitable polynomial, fit using an appropriate re-gression method. It then could be possible to determine a local CBL heightat the point of maximum gradient on this fitted curve.6.2.2 Examine Resolution EffectsRuns could be carried out at lower resolution, to examine the effects on thecurves in Figure 4.19 and eliminate or confirm this as a cause for disagree-ment with the results of Federovich et al. (2004) discussed in Section 5.3.1.The height and θ jump definitions of Garcia and Mellado (2014) could be82matched to facilitate direct comparison with their results and speak moreto the need for increased resolution.6.2.3 Further Explore Entrainment RegimesAs discussed already in this Section, there is sufficient cause to assumethere is a change in entrainment mechanism as Ri increases. Animated vi-sualizations of two-dimensional horizontal slices showed thermals regularlyimpinging on the FA with associated, intermittent periods of vigorous ac-tivity. A possible way to observe these mixing events and how they changewith respect to Ri and time, is to measure turbulent velocity and vorticityat local points, or sub-domain regions. The connection between increasedhorizontal and downward motions in the EZ, and CBL growth can, easily betested by concurrently measuring the local height. Furthermore, turbulentactivity measured at different levels within the EZ could shed further lighton the turbulence characteristics of Garcia and Mellado’s (2014) suggestedtwo layer structure.6.2.4 Apply a Mass Flux SchemeThe robustness of equations 1.16 and 1.25 could be tested by first estab-lishing criteria for identifying CBL air and then calculating the entrainmentrate based on the increase in its volume using the method described in Daweand Austin (2010). CBL air could be identified using a passive tracer, andor potential temperature. The EZ could, for example, be defined based onstatistics of local CBL heights.83ReferencesBatchvarova, E. and S.-E. Gryning, 1994: An applied model for theheight of the daytime mixed layer and the entrainment zone.Bundary-Layer Meteorology, 71, 311–323.Behnel, S., R. Bradshaw, C. Citro, L. Dalcin, D. S. Seljebotn, andK. Smith, 2011: Cython: The best of both worlds. Computing in ScienceEngineering, 13, 31 – 39, doi:10.1109/MCSE.2010.118.Betts, A. K., 1974: Reply to comment on the paper: ”non-precipitatingcumulous convection and its parametrization”. Quart. J. Roy. Meteor.Soc., 100, 469 – 471.Boers, R., 1989: A parametrization of the depth of the entrainment zone.Journal of Applied Meteorology, 107–111.Brooks, I. M. and A. M. Fowler, 2012: An evaluation of boundary-layerdepth, inversion and entrainment parameters by large-eddy simulation.Bundary-Layer Meteorology, 142, 245–263.Crum, T. D. and R. B. Stull, 1987: Field measurements of the amount ofsurface layer air versus height in the entrainment zone. Journal of theAtmospheric Sciences, 44, 2743 – 2753.Crum, T. D., R. B. Stull, and E. W. Eloranta, 1987: Coincident lidar andaircraft observations of entrainment into thermals and mixed layers.Journal of Climate and Applied Meteorology, 26, 774–788.Dawe, J. T. and P. H. Austin, 2010: Interpolation of les cloud surfaces foruse in direct calculations of entrainment and detrainment. MonthlyWeather Review, 139, 444–456.Deardorff, J. W., 1970: Convective velocity and temperature scales forthe unstable planetary boundary layer and for rayleigh convection.Journal of the Atmospheric Sciences, 27, 1211 – 1213.84Deardorff, J. W., 1972: Numerical investigation of neutral and unstableplanetary boundary layers. Journal of the Atmospheric Sciences, 29, 91 –115.Deardorff, J. W., 1979: Prediction of convective mixed-layer entrainmentfor realistic capping inversion structure. Journal of the AtmosphericSciences, 36, 424–436.Deardorff, J. W., G. E. Willis, and B. J. Stockton, 1980: Laboratorystudies of the entrainment zone of a convectively mixed layer. J. FluidMech., 100, 41–64.Ebert, E. E., U. Schumann, and R. B. Stull, 1989: Nonlocal turbulentmixing in the convective boundary layer evaluated from large-eddysimulation. Journal of the Atmospheric Sciences, 46, 2178 – 2207.Federovich, E., R. Conzemus, and D. Mironov, 2004: Convectiveentrainmnent into a shear-free, linearly stratifies atmosphere: Bulkmodels reevaluated through large eddie simulation. Journal of theAtmospheric Sciences, 61, 281 – 295.Fedorovich, E. and R. Conzemius, 2001: Large Eddy Simulation ofConvective Entrainment in Linearly and Discretely Stratified Fluids. 1sted., Kluwer Academic Publishers.Garcia, J. R. and J. P. Mellado, 2014: The two-layer structure of theentrainment zone in the convective boundary layer. Journal of theAtmospheric Sciences, doi:10.1175/JAS-D-130148.1.Hunter, J. D., 2007: Matplotlib: A 2d graphics environment. ComputingIn Science & Engineering, 9 (3), 90–95.Jones, E., T. Oliphant, P. Peterson, et al., 2001–: SciPy: Open sourcescientific tools for Python. URL http://www.scipy.org/, [Online; accessed2014-12-22].Khairoutdinov, M. F. and D. A. Randall, 2003: Cloud resolving model ofthe arm summer 1997 iop: Model formulation results, uncertainties andsensitivities. Journal of the Atmospheric Sciences, 60, 607–623.Kolmogorov, A. N., 1962: A refinement of previous hypotheses concerningthe local structure of turbulence in a viscous incompressible fluid at hightreynolds number. Journal of Fluid Mechanics, 13, 82–85.85Mahrt, L. and J. Paumier, 1984: Heat transport in the atmosphericboundary layer. Journal of the Atmospheric Sciences, 41, 3061–3075.Nelson, E., R. Stull, and E. Eloranta, 1989: A prognostic relationship forentrainment zone thickness. Journal of Applied Meteorology, 28, 885–901.Schmidt, H. and U. Schumann, 1989: Coherent structure of theconvective boundary layer derived from larg-eddy simulations. J. Fluid.Mech., 200, 511–562.Sorbjan, Z., 1996: Effects caused by varying the strength of the cappinginversion based on a large eddy simulation model of the shear-freeconvective boundary layer. Journal of the Atmospheric Sciences, 53,2015–2024.Steyn, D. G., M. Baldi, and R. M. Hoff, 1999: The detection of mixedlayer depth and entrainment zone thickness from lidar backscatterprofiles. Journal of Atmospheric and Oceanic Technology, 16, 953–959.Stull, R. B., 1973: Inversion rise model based on penetrative convection.Journal of the Atmospheric Sciences, 30, 1092–1099.Stull, R. B., 1976a: The energetics of entrainment across a densityinterface. Journal of the Atmospheric Sciences, 33, 1260 – 1267.Stull, R. B., 1976b: Mixed layer depth model based on turbulentenergetics. Journal of the Atmospheric Sciences, 33, 1268 – 1278.Stull, R. B., 1988: An Introduction to Boundary Layer Meteorology. 1sted., Kluwer Academic Publishers.Sullivan, P. P., C.-H. Moeng, B. Stevens, D. H. Lenschow, and S. D.Mayor, 1998: Structure of the entrainment zone capping the convectiveatmospheric boundary layer. Journal of the Atmospheric Sciences, 55,3042–3063, doi:10.1007/s10546-011-9668-3.Sullivan, P. P. and E. G. Patton, 2011: The effect of mesh resolution onconvective boundary layer statistics and structures generated by largeeddie simulation. Journal of the Atmospheric Sciences, 58, 2395–2415,doi:10.1175/JAS-D-10-05010.1.Tennekes, H., 1973: A model for the dynamics of the inversion above aconvective boundary layer. Journal of the Atmospheric Sciences, 30,558–566.86Traumner, K., C. Kottmeier, U. Corsmeier, and A. Wieser, 2011:Convective boundary-layer entrainment: Short review and progress usingdoppler lidar. Bundary-Layer Meteorology, 141, 369–391,doi:10.1007/s10546-011-9657-6.Turner, J. S., 1986: Turbulent entrainment: the development of theentrainment assumption and its application to geophysical flows. J. FluidMech., 173, 431–471.Vieth, E., 2011: Fitting piecewise linear regression functions to biologicalresponses. Journal of Applied Physiology, 67, 390–396.Wilde, N. P., R. A. Stull, and E. W. Eloranta, 1985: The lcl zone andcumulous onset. Journal of Cimate and Applied Meteorology, 24, 640 –657.87A. AppendicesA.1 Potential Temperature: θθ = T(p0p)Rdcp(A.1)p0 and p are a reference pressure and pressure respectively.cpθdθdt =cpTdTdt −Rdpdpdt (A.2)If changes in pressure are negligible compared to overall pressure, as in thecase of that part atmosphere that extends from the surface to 2 km aboveit.cpdθθ = cpdTT −Rdpdpp (A.3)dθθ =dTT (A.4)and ifθT ≈ 1 (A.5)then small changes in temperature are approximated by small changes inpotential temperaturedθ ≈ dT or θ′ ≈ T ′ (A.6)and at constant pressure change in enthalpy (H) is88dH = cpdT. (A.7)This serves as justification for defining w′θ′ as the vertical heat flux.A.2 Second Law of Thermodynamicsdsdt ≥qT (A.8)For a reversible processdsdt =qT (A.9)Using the first law and the equation of state for an ideal gasqT =1T(dhdt − αdpdt)=cpTdTdt −Rdpdpdt (A.10)sodsdt =qT =cpθdθdt (A.11)For a dry adiabatic atmospheredsdt =cpθdθdt = 0 (A.12)A.3 Reynolds Decomposition and Simplificationof Conservation of Enthalpy (or Entropy) fora Dry Atmosphere∂θ∂t + ui∂θ∂xi= νθ∂2θ∂x2i− 1cp∂Q∗∂xi(A.13)ν and Q∗ are the thermal diffusivity and net radiation respectively. If weignore these two effects then∂θ∂t + ui∂θ∂xi= 0 (A.14)89θ = θ + θ′ , θ = ui + u′i (A.15)∂θ∂t +∂θ′∂t + ui∂θ∂xi+ u′i∂θ∂xi+ ui∂θ′∂xi+ u′i∂θ′∂xi= 0 (A.16)Applying Reynolds averaging gives∂θ∂t + ui∂θ∂xi+ u′i∂θ′∂xi= 0 (A.17)Ignoring mean winds∂θ∂t + u′i∂θ′∂xi= 0 (A.18)using flux form∂θ∂t +∂(u′iθ′)∂xi− θ′ ∂u′i∂xi= 0 (A.19)under the bousinesq assumption ∆ · ui = 0∂θ∂t = −∂(u′iθ′)∂z (A.20)ignoring horizontal fluxes∂θ∂t = −∂(w′θ′)∂z (A.21)A.4 Tri-Linear Fit for Determining Local MLHeight hl0The following is a modified version of the piecewise linear regression methodused in Vieth (2011) and was implemented using Cython (Behnel et al.2011). Potential temperature is assumed to be linear function of heightθ = bz + a (A.22)90.Each local θ profile was assumed to have three linear portions, with slopes(b1, b2, b3) and intercepts (a1, a2, a3) as follows:b1 =∑j0 z(i)θ(i)− 1j∑j0 z(i)∑j0 θ∑j0 z(i)2 − 1j (∑j0 z(i))2(A.23)a1 =∑j0 z(i)θ(i)∑j0 z(i)− b1∑j0 z(i)2∑j0 z(i)(A.24)b2 =∑kj z(i)θ(i)− (k − j)a1 + b1z(j)∑kj z(i)− (k − j)z(j)(A.25)a2 =∑kj z(i)θ(i)∑kj z(i)− b2∑kj z(i)2∑kj z(i)(A.26)b3 =∑nk z(i)θ(i)− (k − j)a1 + b1z(j)∑kj z(i)− (k − j)z(j)(A.27)a3 =∑nk z(i)θ(i)∑nk z(i)− b3∑nk z(i)2∑nj z(i)(A.28)where z(i) and θ(i) are a local height and potential temperature value ata particular height index i. j is the height index of the ML top, h0. k isthe height index for the top of the EZ, h1. n is the total number of heightlevels. The best fit is that with the smallest residual sum of squaresRSS(j, k) =j∑0(θ(i)−(a1+b1z(i)))2+k∑j(θ(i)−(a2+b2z(i)))2+n∑k(θ(i)−(a3+b3z(i)))2(A.29).91
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- The effects of two parameters on convective boundary...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
The effects of two parameters on convective boundary layer entrainment Chaparro, Niamh 2014
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | The effects of two parameters on convective boundary layer entrainment |
Creator |
Chaparro, Niamh |
Publisher | University of British Columbia |
Date Issued | 2014 |
Description | The atmospheric convective boundary layer has been studied for over thirty years in order to understand the dynamics and scaling behaviour of its growth by entrainment. This enables prediction of its entrainment rate and entrainment zone depth, and so parameterizations thereof for use in global circulation models. Fundamentals, such as the dependence of the entrainment rate and entrainment zone depth on the convective Richardson number, have been established but there is still unresolved discussion about the form of these relationships. Details regarding the structure of the entrainment zone continue to emerge. The variety of convective boundary layer height and entrainment zone depth definitions adds further complexity. The study described in this thesis aims to join this ongoing discussion. A dry, shear-free, idealized convective boundary layer in the absence of large scale winds was modeled using a large eddy simulation. The use of ten ensemble cases enabled calculation of true ensemble averages and potential temperature fluctuations as well as providing smooth average profiles. A range of convective Richardson numbers was achieved by varying the two principle external parameters: surface vertical heat flux and stable upper lapse rate. The gradient method for determining local convective boundary layer height was found to be unreliable so a multi-linear regression method was used instead. Distributions of the local heights thus determined were found to narrow with increased upper stability. Height and entrainment zone depth were then defined based on the ensemble and horizontally averaged potential temperature profile. The resulting relationships of entrainment rate and entrainment zone depth to Richardson number showed behaviour in general agreement with theory and the results of other studies. The potential temperature gradient in the upper convective boundary layer and entrainment zone was seen to depend on the upper lapse rate, as was the positive downward moving temperature fluctuations at the \acs{CBL} top. Overall, once the surface heat flux was accounted for by applying the \acs{CBL} height as a scale, the upper lapse rate emerged as the dominant parameter influencing scaled entrainment zone depth, and potential temperature variance in the entrainment zone and upper convective boundary layer. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2015-01-06 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivs 2.5 Canada |
DOI | 10.14288/1.0167268 |
URI | http://hdl.handle.net/2429/51761 |
Degree |
Master of Science - MSc |
Program |
Atmospheric Science |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2015-02 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/2.5/ca/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2015_february_chaparro_niamh.pdf [ 5.68MB ]
- Metadata
- JSON: 24-1.0167268.json
- JSON-LD: 24-1.0167268-ld.json
- RDF/XML (Pretty): 24-1.0167268-rdf.xml
- RDF/JSON: 24-1.0167268-rdf.json
- Turtle: 24-1.0167268-turtle.txt
- N-Triples: 24-1.0167268-rdf-ntriples.txt
- Original Record: 24-1.0167268-source.json
- Full Text
- 24-1.0167268-fulltext.txt
- Citation
- 24-1.0167268.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}]}"
data-media="{[{embed.selectedMedia}]}"
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.24.1-0167268/manifest