Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Numerical simulation of canopy flow and carbon dioxide flux at the west coast flux station Sun, Haizhen 2005

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

Item Metadata

Download

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

Full Text

N U M E R I C A L S I M U L A T I O N OF C A N O P Y F L O W A N D C A R B O N DIOXIDE F L U X AT T H E W E S T C O A S T F L U X STATION  by  HAIZHEN SUN B.Sc, Nanjing University, 1992  A THESIS S U B M I T T E D IN P A R T I A L F U L F I L M E N T OF T H E R E Q U I R E M E N T S F O R T H E D E G R E E OF M A S T E R OF SCIENCE in  T H E F A C U L T Y OF G R A D U A T E STUDIES  (Atmospheric Science)  T H E U N I V E R S I T Y OF BRITISH C O L U M B I A  April 2005  © Haizhen Sun, 2005  Abstract Using an eddy covariance-approach to estimate net-ecosystem-exchange (NEE) of carbon dioxide (CO2) requires that advection is negligible. This approach assumes flat and homogenous terrain. Sloping and heterogeneous topography may contribute to errors in estimating carbon sequestration or loss of an ecosystem. The complex topography at the West Coast Flux station on Vancouver Island raises uncertainties on the estimation of NEE. This research employs a high-resolution atmospheric model to study the effects of a deep forest canopy on atmospheric boundary-layer flow, and to evaluate the roles of mesoscale flow and the advection effects on NEE estimates. The model is refined to include tree-drag, radiation effects of canopy on surface energy budget, C 0  2  budget and soil conduction. Over an idealized  two-dimensional mountain of similar horizontal and vertical scales as Vancouver Island, the numerical simulations capture the first-order effects of diurnal heating/cooling on the sloping terrain during fair-weather. The volume exchange of heat, momentum and CO2 above and within the canopy appears to be strongly affected by the local flows resulting from diurnal thermal forcing, land/sea breezes and convective thermals. Synoptic-scale forcing is neglected. The simulated C 0 concentration shows significant variation and gradients near the surface, 2  corresponding to the diurnal change of stability in the atmospheric boundary layer. The time evolution of C 0 spatial field exhibits the transport processes at night by drainage flows. The C 0 2  2  budget analysis further quantifies the contributions of horizontal and vertical advection. The estimated NEE by only using the sum of the storage change and vertical flux results in large fluctuation during daytime and underestimation of nocturnal respiration by about 40%. When we reduce the effects of strong convection by subtracting a half-hour average from instantaneous variables, as is done in eddy-covariance approach, the daytime fluctuations of estimated N E E are strongly decreased and the estimated nocturnal respiration is within 20% of the prescribed source. Back-streakline analysis shows that the advection source areas of C 0 reaching the flux towers 2  during fall are mainly from the northwest and southeast along Georgia Strait with a larger advective source area during daytime and a smaller one at night.  Table of Contents Abstract  ii  Table of Contents  iii  List of Tables  v  List of Figures  vi  List of Symbols and Abbreviations  ix  Acknowledgements  xii  Chapter 1 Introduction and research goals  1  1.1 Background 1.2 The West Coast Flux Station 1.3 Eddy covariance method and NEE estimates 1.4 Questions of heterogeneity and non-stationarity 1.5 Goals of this research Chapter 2  Refinement of the Weather Fire Integrated System (WFIS)  2.1 Introduction 2.2 Treatment of canopy drag 2.2.1 Basic variables describing canopy 2.2.2 Treatment of canopy drag 2.3 Refinement of model physics 2.3.1 Energy budget within a canopy 2.3.2 Soil surface heat budget 2.4 Summary Chapter 3  Simulations of canopy flow over a hill  3.1 Introduction 3.2 Experimental design 3.3 Vertical grid design 3.4 Results and discussion 3.4.1 Surface temperature variations 3.4.2 Katabatic and anabatic flow 3.4.3 Drainage flow 3.4.4 Thermals and land/sea breezes 3.4.5 Vertical profiles 3.5 Conclusions Chapter 4  1 2 3 7 9  C 0 budget and advection effects 2  4.1 Introduction 4.2 Treatment of C 0 budget in the model 4.3 C 0 budget experiments and results 2  2  11  :  11 13 13 15 17 17 20 21 23 23 26 28 30 30 32 35 36 40 43 46 46 48 52  4.3.1 Simulated soil temperature 4.3.2 Simulated C 0 flux by respiration and photosynthesis 4.3.3 Simulated C 0 profiles in the A B L 4.3.4 Diurnal variations of C 0 concentration and gradients 4.3.5 C 0 spatial distribution 4.4 C 0 budget analysis 4.4.1 Contributions of each term in the C 0 budget 4.4.2 Contributions of each term in NEE estimation 4.5 Summary and Conclusions 2  2  2  2  2  2  Chapter 5 5.1 5.2 5.3 5.4 Chapter 6  Exploratory analysis of the advection source area Introduction Methodology Results Conclusion Summary, Conclusions, Limitations and Future work  6.1 Summary and Conclusions 6.2 Limitations 6.3 Future work Bibliography  80 80 81 84 85 89 89 93 94 96  Appendix A WFIS model dynamics A . l Momentum equations in spherical geometry A.2 Continuity equation A. 3 Terrain-following coordinates Appendix B WFIS model physics B. l Conservation equation for potential temperature B. 2 Radiation cooling schemes Appendix C  53 54 55 58 63 65 66 71 77  Sensitivity experiments and preliminary results  C. l Sensitivity to tree-drag coefficient C.2 Sensitivity to topography and upstream flow  105 105 105 105 107 107 107 110 110 115  List of Tables 3.1  Outlines the characteristics of diurnal-cycle experiments  28  4.1  Characteristics of C 0 budget experiments  52  c.l  Characteristics of sensitivity to drag-coefficient experiments  Ill  c.2  Characteristics of sensitivity to topography and to upstream flow experiments  116  2  List of Figures 1.1  The global average atmospheric carbon dioxide mixing ratio and the long-term trend of growth between 1000 and 2000  1  1.2  The Mature Stand site at the West Coast Flux Station  3  1.3  The time series of measurements of wind speed, wind direction, net radiation, and air temperature (d), CO2 eddy flux (e) at main tower site of West Coast Flux station on Vancouver Island during days 207-214, 2003  8  2.1  Vertical distribution of leaf area density  16  3.1  Schematic diagram of topography and canopy specifications for the outer domain (upper) and the nested domain (lower)  27  3.2  Vertical grid structure for the two nests  29  3.3  Temperature of soil superficial layer (2.5 cm, thick line) and surface air temperature (3 m, thin line) for NOF0 {r\ = 0), NOF1 (r\ = 0.8), and NOF2 (r) = 1)  31  3.4  Similar to Fig. 3.3 but for experiments INF [0-2] with dynamics on  32  3.5  Horizontal wind at the four times of 0600, 1200, 1800 and 0000 LST in the outer nest for INF2. Contour interval is 1 m/s with solid showing positive values from west to east and dashed showing negative  33  3.6  Similar to Fig. 3.5 except for the four times of 1730, 1745, 1800 and 1830 LST  34  3.7  Drainage flow for experiments INFO and INF2  35  3.8  Thermal structure and land-sea breezes at 1400 and 0200 LST for nested domain for experiments INF2  37  3.9  Wind field at 1400 LST for the lower 150 m of Fig. 3.8  .....38  3.10  Statistics of horizontal (a) and vertical (b) velocity at 6 m above ground, and horizontal (c) and vertical (d) at 42 m above the ground at grid 280 for experiments INF2 (r\ = 1.0) 39  3.11  Potential temperature profiles for daytime (left) and nighttime (right). The lower panels show air potential temperature from the ground to 50 m  41  3.12  Similar to Fig. 3.11 except for horizontal wind  42  3.13  Similar to Fig. 3.12 except for vertical velocity  43  4.1  Schematic diagram of carbon sources/sinks in an ecosystem  50  4.2  Radiation budget at the canopy top for a 4-day simulation NOF2  51  4.3  Variation of simulated soil temperatures with time at several soil depths for r\ = 1.0 case and r| = 0 case  53  4.4  Simulated soil respiration for experiment C02F1 and CO2F0  54  4.5  Simulated photosynthetic C 0  2  flux using empirical formula (4-6) (dash line) and  parameterized scheme (4-7) (solid line) 4.6  54  Simulated vertical profiles of C 0 concentration at 6 times (LST) during a day at grid 2  location 392  56  4.7  Similar to Fig. 4.6 except for grid point 450  57  4.8  The time series of C 0 concentration (upper panel) at the height of 6 m, 18 m and 33 m, and 2  C 0 concentration gradient between the three layers (lower panel) for the outer domain of the 2  dense-canopy experiment C02F1 4.9  60  Diurnal courses of C 0 concentration at the height of 3 m, 9 m and 15 m, and C 0 gradient 2  2  between these three vertical layers for nested domain of the dense-canopy experiment C02F1  61  4.10 Same as Fig.4.8 except for the bare-ground case (CO2F0)  62  4.11 Time series of horizontal C 0 gradient at 3 m high above ground at grid point 280 for the 2  dense-canopy case  62  4.12 Simulated C 0 spatial distributions at different times in a day for dense-canopy experiment 2  C02F1  64  4.13 Simulated contribution of each term in the C 0 budget equation (4-10) at 1300, 1900 and 2  0100 LST  ;  68  4.14 The contributions of horizontal advection (squares, term II+IV), vertical advection (asterisks, term III+V) and storage (circles, term I) at 1300, 1900, and 0100 LST  69  4.15 Estimation of the net source/sink using the sum of different terms  70  4.16 Diurnal cycle of each term in the NEE estimation (4-12)  74  4.17 Time series of the contributions of horizontal and vertical advections  75  4.18 Comparison of the incorporated source to the estimated source using: a) the sum of five terms and b) the sum of three terms  76  4.19 Comparison of the incorporated source to the estimated source using a similar analysis procedure with the EC approach  76  5.1  MC2 model domain configuration  '...  ..81  5.2  Surface wind and back streaks with trajectories passing through the northern flux site at 16  pm, Dec. 29, 2004 5.3  The locations of air parcels 8 hours before passing through the Northern site at 1600 LST in October (a), November (b) and December (c), and Oct. to Dec. (d) 2004  5.4  83  86  The locations of air parcels 8 hours before passing through the Northern site at 4 am PST in October (a), November (b) and December (c), and Oct. to Dec. (d) 2004  87  5.5  Same as Fig.5.3 except for the Southern (Young Fir) site  87  5.6  Same as Fig.5.4 except for the Southern (Young Fir) site  88  C. 1 Streamlines for experiments DRAG0D and DRAG ID for the outer domain  112  C.2  Inner domain streamlines for experiments DRAG0D (a) and DRAG1D (b)...'.  112  C.3  Horizontal velocity for inner nest for experiment DRAG1 [A-D]  114  C.4  Vertical velocity for inner nest for experiment DRAG 1 [A-D]  115  C.5  Horizontal velocity for the third nested domain for experiments HMNT [1-4]  118  C. 6 Horizontal velocity for the second nest for RGMV [ 1 -4] experiments  118  Symbols and Abbreviations Greeks air density canopy density quantum yield ground albedo canopy albedo emissivity of the canopy emissivity of the ground surface stress a fraction of area covered with trees air potential temperature latitude longitude declination of the Sun extinction coefficient finite-difference operator horizontal grid size vertical grid size  p p a ct a e e x a  c  G  t  c  G  11 4> ^ 5  r ^ ^  z  ^ ^-gL . s A  A aerodynamic surface area of plant atmospheric boundary layer leaf area density canopy photosynthetic capacity B  B  Bowen Ratio  co c c' Q  2  Com  £ J  M m  1  C carbon dioxide mixing ratio of carbon dioxide fluctuation of carbon dioxide concentration surface drag coefficient drag coefficient dry mass of the canopy per unit ground area air specific heat moist mass of the canopy per unit ground area canopy specific heat -l 1 , • 1 soil volumetric heat capacity  D Def  D molecular diffusion deformation  g£  E eddy covariance  F drag force eddy covariance flux soil efflux vertical C 0 flux  F F F F  c 0 v  2  G G  G heat conduction from the ground fraction of solar radiation  h H Hq HR  H canopy height sensible heat flux within the canopy sensible heat flux to the ground solar hour angle  K k k K K  K eddy-mixing coefficient soil thermal diffusivity extinction coefficient heat transfer coefficient water vapor transfer coefficient  s  m  s  H w  LAI L(z) LE LEg L  L cumulative leaf area index local leaf area index latent heat flux in the air latent heat flux from the soil latent heat of condensation or vaporization  N NEE NWP  N Brunt-Vaisaila frequency net ecosystem exchange numerical weather prediction  P P P ppm  P surface pressure energy consumed via photosynthesis carbon dioxide consumed by photosynthesis parts per million  Q (2i  Q Photosynthetic active radiation Increase factor of R for a 10 °C increase of T  Y  RS S  0  s  s  Ri i?NH 7?Lh4i? l R i? R i? R  R Local Richardson number net radiation at the treetop incoming long-wave radiations outgoing long-wave radiations net radiation within the canopy respiration rate at a reference temperature T soil respiration rate net radiation without canopy net radiation to the ground  SGS S S  S Sub-grid-scale direct solar radiation incoming radiation at the top of the atmosphere  j 7 7 X  T air temperature a reference temperature soil temperature ground surface temperature  Lh  Np  ref  s  N  NC  0  ref s  0  u  <  u  v  i , -p.. y  y  ref  U wind velocity in x direction wind fluctuation in x direction V wind velocity in y direction wind fluctuation my direction wind vector wind velocity at the top of canopy  w  <  W wind velocity in z direction wind fluctuation in z direction  x  X Cartesian coordinate in north-south direction  y  Y Cartesian coordinate in west-east direction  w  Z z  Z zenith angle of the sun Cartesian coordinate in the vertical direction  Acknowledgements I would first like to express my gratitude to my supervisor, Dr. R.B. Stull for his continuous support, encouragement and intellectual guidance. Special thanks goes to committee member Dr. T.L. Clark for his invaluable help and endless streams of ideas. I must also thank committee member, Dr. T.A. Black who is always available for timely advice and insight. It is their precious support and help that led to the completion of this thesis. I gratefully acknowledge the advice and observation data provided by Dr. Morgenstern, Dr. Jassal and Tiebo Cai in the Biometeorology/Soil physics Group at UBC. I would also like to recognize Dr. Martilli and Dr. Modzelewski for their valuable comments and recommendation on literatures. My thanks go to all the people in the Weather Forecast Research Team that have made the 2+ years here so memorable: enjoyable lunch-times and discussion with Yongmei Zhou and Xingxiu Deng; pleasurable English-learning times with Trina Cannon; and delightful discussion with George Hicks II, Luca Delle Monache and Yan Shen. Last, but not least, my greatest debt of gratitude goes to my husband Aiming, my daughter Zina and my mother for their unwavering support and understanding. This work is supported by Fluxnet Canada Research Network (FCRN), based on grants from the Canadian Foundation for Climate and Atmospheric Science (CFCAS), and the Natural Sciences and Engineering Research Council (NSERC) of Canada. Additional support for computer equipment and support personnel came from Environment Canada, the BC Ministry of Water, Land and Air protection, and the BC Forest Innovation Investment account. Computers in the U B C Geophysical Disaster Computational Fluid Dynamics Center used for this study were purchased with grants from the Canadian Foundation for Innovation, the BC Knowledge Development Fund, and UBC. System administration for these computers was supported by the NSERC C3 Technical Analyst support program, and with additional infrastructure funds from NSERC via the U B C Faculty of Science.  xu  Chapter 1 Introduction and research goals 1.1 Background Carbon dioxide (C0 ) 2  is an important greenhouse gas affecting the Earth's climate.  Atmospheric CO2, which has increased since the beginning of the industrial revolution from around 280 ppm (parts per million) in 1800 to almost 370 ppm today (Fig. 1.1), is significantly influenced by human activities [Schimel, 1995]. This C 0 increase enhances 2  the greenhouse effects of the atmosphere, and likely contributes to global warming and climate change [Cox, 2000]. Therefore understanding and quantifying the global carbon budget is of considerable importance i f the magnitude of possible future global warming is to be accurately predicted.  '  YEAR  Figure 1.1: The global average atmospheric carbon dioxide mixing ratio and the long-term trend of growth between 1000 and 2000. [After Sarmiento and Gruber, 2002] The global transport of carbon (partly in the form of CO2) among the land, oceans and the atmosphere is called the Global Carbon Cycle. The CO2 emitted into the atmosphere, together with the uptake by the terrestrial sinks and oceans, governs the CO2 content of the atmosphere, which is observed by the global sampling network. Estimates of the total oceanic uptake by different methods seem to converge [Fan et al., 1998], and the remainder  must be absorbed by "sink" processes on the land. But we lack a quantitative understanding of where the sinks are, how they work, how long they will keep working, and how they can be altered by human activities. Therefore, direct and long-term carbon flux measurements and model calculations are needed to clarify the roles of sources or sinks of different regions situated under different climatic conditions, and to investigate the behavior of different ecological systems. In response to the question of the 'missing carbon sink', interest in forest microclimate and exchange mechanisms between forests and the atmosphere has blossomed over the last decade. Continental networks such as EuroFlux [Valentini et al., 2000], AmeriFlux [Baldocchi et al., 2001] and global networks such as Fluxnet [Baldocchi et al., 2001] have been built to monitor annual carbon uptake, water vapor and energy between terrestrial ecosystems and the atmosphere [http://www-eosdis.ornl.gOv/FLUXNET/1. The goal of Fluxnet is to understand the mechanisms controlling the exchanges of CO2, water vapor and energy between ecosystems and the atmosphere across a spectrum of time and space scales.  1.2 The West Coast Flux Station At present, over 200 flux tower sites are operating on a long-term and continuous basis worldwide. The West Coast Flux Station, also called 'British Columbia Flux Station', is one of seven flux stations of Fluxnet-Canada rhttp:/www.fluxnet-canada.ca]. It includes three tower sites: Mature Stand, Young Fir and Clear-cut. The Mature Stand site (shown in Fig. 1.2) is located on the east side of Vancouver Island near Campbell River on a 5-10 degree slope oriented east/northeast. Georgia Strait, an inland waterway, is 9 kilometers east of the site. The site is covered with a typical second-growth forest. The Clearcut site. (49.958° N , 125.433° W) is close to the Mature Stand site. The Young Fir site (49.519° N , 124.902° W) is much closer to the Georgia Strait and is covered with 13-14 year old Douglas Fir. As most of our numerical experiments are performed over a domain configured to mimic the west-to-east cross-section of Vancouver Island through the Mature Stand site, here we list the characteristics of the Mature Stand site:  •  Latitude: 49.905° N  •  Longitude: 125.366° W  •  Climate: temperate  •  Elevation: 320 m  •  History: -56 years old Douglas fir  •  Canopy height: 33 m  •  Measurements: continuous  •  Tower height: 45 m  •  Flux measurement height: 42 m  •  Vegetation  type:  Temperate  Douglas fir with 17% red cedar and 3% western hemlock  Figure 1.2: The Mature Stand site at the West Coast Flux Station [from http://www.fluxnetcanada.ca]  Since 1997 the U B C Biometeorology/Soil Physics Group has been monitoring radiation, heat, moisture, and carbon dioxide fluxes at the West Coast Flux station using an eddy covariance technique (EC) applied above forest canopy. The E C measurements of atmosphere/biosphere vertical exchange of CO2 have been used to estimate the net ecosystem-atmosphere exchange of CO2 (NEE) to further examine biophysical processes that control carbon release or uptake by the vegetation [Drewitt et al., 2002; Humphreys et a l , 2003; Morgenstern et al., 2004].  1.3 Eddy covariance method and NEE estimates The eddy covariance (EC) method, a technique to measure turbulent fluxes, was first used about 40 years ago [Baldocchi, 2003]. But it is only recently that technological advances have permitted the E C method to be used to assess ecosystem-atmosphere exchange of heat, momentum, moisture and chemical compounds such as carbon dioxide. Compared to the chamber technique to estimate N E E , the E C technique has two desirable properties: it does not disturb the vegetation and it can be used for continuous long-term operation. Both of these are clear advantages over chambers [Davidson et al., 2002]. Furthermore, it produces  a direct measure of net CO2 exchange across the canopy-atmosphere interface [Tans et al., 1996; Baldocchi, 1988; 2003]. The atmosphere contains turbulence motions of upward and downward moving air that transport scalars (such as CO2, H2O and heat) and momentum. The E C technique samples these turbulent motions using fast-response sensors to determine the net vertical transport of the above entities moving across the canopy-atmosphere interface. Practically, this task is accomplished by statistical analysis of the instantaneous vertical flux measurements using Reynolds' rules of averaging [Webb et al., 1980; Stull, 1988]. The result of the operation is a relationship that expresses the mean flux density (averaged over some time span such as half-hour or one-hour) as the covariance between the vertical velocity (w) and the entity. For CO2, this relationship is: F = p W? c  (1-1)  a  where c is the mixing ratio of C 0 , p is the density of dry air, the over-bars indicate the 2  a  mean value for the time span, and the primes indicate the fluctuations from the mean (i.e., c' = c-c).  A l l symbols are defined on page ix to xi.  The C O 2 conservation equation provides the basic theoretical framework  for  implementing the eddy-covariance technique to deduce the source/sink strength of an ecosystem [Baldocchi et al, 1988; Massman et a l , 2002; Baldocchi, 2003] and to interpret micrometeorological flux measurements. This equation is: ct  dt  cy  ,ck  dc  cy  &  c  where S is the sum of C 0 sources and sinks in a control volume, and u, v, and w are the c  2  wind velocity components along x, y and z directions, respectively. D represents the molecular diffusion term, which can be neglected since it is several orders of magnitude smaller than the turbulence diffusion term in the atmosphere. By ignoring molecular diffusion, (1-2) can be written as: d (I)  &  d)> (II)  &  dc  cy (III)  ck (IV)  c (V)  Equation (1-3) states that (I) the time rate of change of CO2 at a fixed point is balanced by mean horizontal and vertical advection (II), by the mean horizontal (III) and vertical (IV)  divergence or convergence of turbulent fluxes and by sources and sinks (V). The C O 2 concentration in a control volume will remain unchanged i f the flux of CO2 entering the volume equals that leaving it. Flux divergence of CO2 acts to remove CO2 from the volume, whereas flux convergence results in the accumulation of CO2 in the controlled volume. If the underlying surface is horizontally homogenous and is on flat terrain (in x and y direction), then terms (II) and (III) equal to zero. The net result is that (1-3) simplifies to a balance between the storage change, the vertical turbulent flux divergence of CO2 and its biological source/sink strength, S (term V): c  m i ^ ^ ' K s  (Z  (1-4)  )  Because above the canopy no biological source/sinks exist, S (z) = 0, and within the c  canopy, usually S {z) * 0, thus, by integrating (1-4) from z = 0 to z = h, where h is canopy c  height, one derives an equality between the mean vertical eddy flux measured at some height above the canopy, p w'c'\ , a  the storage change, and the net flux of material in and  h  out of the underlying soil F and vegetation, and the integrated source effects in the canopy 0  space: Po^'  l„  +Pa )%dz = F +) 0  0  S ( )dz c  Z  (1-5)  0  0  h  In practice, it is the terms, p w'c'\ a  h  and p j(8c/dt)dz, that are evaluated in using the E C a  o  technique. Equation (1-5) forms the basic theoretical framework for using the E C method to estimate sources/sinks of scalars [Baldocchi, 2003]. The E C technique has been widely used in observational studies of CO2 sources/sinks for different terrestrial ecosystems [Wofsy et al., 1993; Black et al., 1996; Saigusa et al., 1998; Holland et al., 1999; Aubinet et al., 2000; Carrara et al., 2003], and of water vapor and sensible heat flux between the biosphere and the atmosphere [Stannard et al., 1994; Anthoni et al., 1999; Twine et a l , 2000], The most popular application of the E C technique is to estimate N E E in response to the current interests in closing the global carbon budget. N E E is defined as the net exchange rate of CO2 between the ecosystem (including the forest and the soil) and the atmosphere in a controlled volume. A s stated by Staebler and Fitzjarrald, [2004], N E E is defined by:  NEE = F +\S dz 0  (1-6)  c  o  where F is the soil efflux of CO2 entering (or leaving) the control volume at the bottom, 0  and the integral describes the sum of all sources and sinks in the canopy space. Positive values of N E E represent upward flux into the atmosphere, or a carbon loss from the ecosystem. Since N E E reflects the change in carbon stored in the ecosystem, it is used to infer whether the ecosystem is a carbon 'sink' or 'source'. Comparing (1-5) and (1-6), one may conclude that N E E can be directly evaluated by the EC  flux  change of C O 2 , i.e., usually  measurements with air column storage h  NEE = p w'c'\ +p \—dz • But one thing we have to keep in mind is that (1-5) is gained by a  h  a  0  assuming that the underlying surface is horizontally homogenous and flat. However, since horizontal homogeneity is unrealistic, one can substitute for  )s dz c  from (1-3) into (1-6) to  give: NEE=F +)s dz 0  c  — = }^t- dz+ L  0  #  J  °  i  dx  ^+\  '  'dz+(pwc),+ J  Id/  —  "  „  1 dz]+(p.W{O  (1-7)  J  dx  „  t  3/  where the integral of the vertical turbulence flux divergence produces a flux at z = 0 which cancels the soil efflux Fo- The vertical velocity (w) is equal to zero at the ground so the integration of mean vertical advection reduces to (p wc) . F o r two dimensional (y, z) flow, 0  h  rearranging (1-7) yields: NEE= ) m b  [I]  +  } ^ p d  [II]  z  +  ^  )  h  [III]  +  ) ^ l d z  [IV]  ]  +  ^  X  (1-8)  [V]  N E E should be calculated using the full equation (1-8). The eddy covariance term [V] is found by using the E C technique, and the first term [I] is significant on short time scales (less than a day), fortunately it can be determined from routine measurements of the CO2 profiles measured from different height under the tower. The remaining terms are much more difficult to measure, but can be estimated from numerical simulation. The goal of this thesis is to use high-resolution numerical modeling to estimate the other terms in (1-8), to yield a better measure of N E E .  1.4 Questions of heterogeneity and non-stationarity The storage change in (1-8) becomes significant on short time scales (less than a day). The second, third and fourth terms in (1-8) can be non-zero because of undulating terrain and patchy vegetation. It is estimated that 1/3 of North America is topographically too complex for eddy-flux techniques to work [Schimel et al., 2002]. Furthermore, most flux tower sites have problems closing the observational energy budget based on the simplified E C method that considers only terms I and V . Direct measurements of the product of vertical velocity fluctuations (w'), temperature fluctuations (T) and the water-vapor mixing ratio fluctuations (r') yield a direct estimate of sensible heat flux (H) and latent heat flux (LE), but the estimated sensible and latent heat fluxes do not add up to the available energy, given by net radiation minus heat storage and heat conduction by soil. Twine et al. [2000] found the E C method underestimates values of H and LE by 10 to 30%, the largest residual to energy budget closure during the Monsoon'90 experiment occurred at the poorest-rated site (in terms of flatness and slope), while the smallest residual was 2 W m" , occurred at the best-rated site [Stannard et al., 1994]. This 2  raises the question of how does non-ideal terrain affect the estimation of N E E ? Observations have reported that the current method used to estimate N E E ( N E E = turbulent flux plus air column storage term) underestimates nocturnal respiration of CO2 on calm nights under low turbulence conditions compared to the chamber measurements [Wofsy et a l , 1993; Black et al., 1996; Lavigne et al., 1997; Mahrt, 1998; Morgenstern et al., 2004]. In most forests, reasonable N E E estimates are obtained during greater intensities of turbulence. However, under conditions of low turbulence, eddy fluxes are significantly smaller although there are no ecological reasons to expect that respiration should be inhibited. Under these conditions, respired CO2 must be transported away from its local source by some other mechanisms undetected by the simplified measurement techniques. From (1-3), potential candidates for these other mechanisms include the mean vertical and horizontal advective flux gradient, and mean horizontal turbulent flux divergence. Therefore there is an urgent need to quantify these atmospheric transporting mechanisms before credible N E E estimates can be gained. But how can we quantify these regional scale mechanisms using a single tower? To measure these mesoscale factors would require a prohibitively expensive array of tall towers, each with E C sensors at many heights.  The E C flux tower at the West Coast Flux station is located on a 5-10 degree slope oriented east/northeast. The Georgia Strait, 9 kilometers east of the site, contributes to a different aerodynamic surface roughness. Therefore the real conditions at the West Coast Flux station are heterogeneous, and the complex topography raises uncertainties on E C measurements and N E E estimates. Is there any advection effect on N E E estimates at the West Coast Flux station? How large are the non-turbulent fluxes of CO2 that are not captured by the traditional E C method? What mesoscale mechanisms mainly contribute to the advection effects? What are the strengths and limitations of tower data analysis procedures currently used to estimate N E E ?  214 Figure 1.3: The time series of measurements of wind speed (a), wind direction (b), net radiation (c), air temperature (d), C 0 eddy flux (e) at main tower site of West Coast Flux station on Vancouver Island during days 207-214, 2003 (from UBC Biometeorology/Soil Physics Group). The arrows in panel (e) indicate the anomalous values of eddy flux around sunset. 2  In addition, the eddy flux measurements of CO2 at the West Coast Flux station show anomalies (a significant peak of upward CO2 flux) just after sunset (arrows in Fig. 1.3e) and the wind direction changes frequently under fair-weather condition (Fig. 1.3b). So what  particular controlling mechanisms cause these regular anomalies appearing in the flux measurements just after sunset at the West Coast Flux station? Answering this question of non-stationarity requires a further understanding of wind and turbulent regimes at the WestCoast Flux station, since wind flow is an important factor for scalar fluxes (heat, water vapor, carbon dioxide, pollutants, etc) within a forested canopy. Therefore a search for a more robust physical understanding of the fluid dynamics and exchange mechanisms is in order at the West Coast Flux station. The difficulties in measuring the mesoscale flows by a single tower suggest that numerical simulation might provide insights into the physics of the mesoscale controlling mechanisms at the West Coast Flux station. Some of the deficiencies of field experiments could be overcome by using high-resolution numerical model to quantify the regional scale wind and CO2 fields within and above a canopy located on complex terrain. In this way, heterogeneous and non-stationary effects can be included to yield a better N E E estimate.  1.5 Goals of this research The first objective of this study is to refine an existing mesoscale atmospheric model, WFIS (Weather Fire Integrated System), so that it properly treats canopy flow and CO2 fluxes at the West Coast flux station. This refinement of the WFIS model includes the following new parameterizations for: (1) the tree drag force induced by canopy elements in the momentum equations; (2) the radiation effects of the canopy on the surface and canopy energy budgets; (3) soil heat conduction; and (4) the CO2 budget with a treatment of the source/sink terms, i.e., canopy photosynthesis and soil respiration. The second goal is to provide insights on the influence of a deep forest canopy on the atmospheric boundary layer ( A B L ) flow and canopy flow over complex terrain in the context of a diurnal cycle, and to understand the mesoscale controlling mechanisms that account for the CO2 flux anomalies (Fig. 1.3e) at the West Coast Flux station. Particularly, we will see how the thermally driven mesoscale circulations affect the flux measurements on sloping terrain. The geography at the West Coast flux station is well suited for this study because of the sloping terrain, the land/sea breeze pattern associated with the shorelines, and the thermally induced mountain and valley flows.  The final objective is to address the underlying advection effects on the E C flux measurements and the estimation of N E E at the West Coast Flux station. By incorporating the CO2 conservation equation into the high-resolution WFIS model, and solving it in the context of the diurnal cycle over a two-dimensional mountain of similar horizontal and vertical scale as Vancouver Island, we will assess the contributions of the horizontal and vertical advection to the N E E estimates. This allows us to numerically assess the strengths and weaknesses of current simplified procedures (NEE = storage term + vertical turbulent flux) used to estimate N E E at complex sites such as the West Coast Flux station. Also, a climatology of advection source areas affecting the flux tower sites will be built by using the back-streakline products from an operational weather prediction model (MC2) to aid in the interpretation of CO2 concentration measurements. The thesis is organized as follows. The refinement of the model is described in Chapter 2. Numerical experiments with the results are discussed in Chapter 3, which focus on the effects of a deep forest canopy on A B L flow and canopy flows in the context of a diurnal cycle over sloping terrain. Chapter 4 uses the resulting refined model to simulate the diurnal response of CO2 field and to find the contributions of each term in the C O 2 budget, allowing N E E to be assessed. In Chapter 5, the advection source areas for Mature Stand and Young Fir sites will be identified for one season from October to December using backward streaklines calculated daily from the surface wind fields generated by a real-time weather prediction model M C 2 . The last chapter summarizes the results, discusses the limitations of this research, and outlines directions for future research.  Chapter 2 Refinement of the Weather Fire Integrated System (WFIS) In this chapter, the basic characteristics of the WFIS model are introduced, followed by the detailed description of the treatment of canopy drag dynamics in WFIS in Section 2.2. The refinement of the model physics, such as the effects of the forest canopy on radiative contributions to the surface and canopy energy budgets, is described in Section 2.3. The last section summarizes the methodology.  2.1 Introduction The numerical model adopted here to incorporate canopy dynamics and thermodynamics is a considerably modified version of the three-dimensional code of Clark [1977]. It is a limited area, anelastic ( V ( p F ) = 0), finite-difference model in which the equations of o  motion and the first law of thermodynamics are solved in a domain that has an irregular lower boundary. The first version of the mesoscale atmospheric model was developed by Clark [1977] during his tenure at the A E S (Atmospheric Environment Service) now called the M S C (Meteorological Service of Canada), and was also known as the Clark model. Cloud physical parameterizations were later added by Clark [1979]. Interactive grid nesting was added by Clark and Farley [1984] and later refined by Clark and Hall [1991, 1996]. With the recent developments of Message Passing Interface (MPI) and fire dynamics, the model has become known as WFIS (Weather Fire Integrated System). WFIS incorporates the non-hydrostatic and anelastic approximations [Clark, 1977; 1979]. Scale analysis indicates that the hydrostatic approximation breaks down when the resolvable horizontal and vertical scales are almost the same order. Thus, non-hydrostatic dynamics are needed to investigate local scale motions. WFIS invokes the 'anelastic' approximation [Ogura and Phillips, 1962; Lipps and Hemler, 1982], in which sound waves are filtered by appropriate modification of the mass continuity equation (see Appendix A ) . To deal with an irregular lower boundary, such as mountains and canopies, WFIS employs a vertically stretched terrain-following coordinate (Appendix A ) . The model can  be configured in either two or three spatial dimensions, and uses two-way interactive grid nesting, which can telescope down to resolve small-scale features in the vicinity of an observation tower. Also the model contains an explicit treatment of microphysics using the Kessler [1969] warm rain and the Koenig-Murray [1976] ice parameterizations. Sub-grid scale effects are closed using the first-order Smagorinsky-Lilly [Smagorinsky 1963, L i l l y 1962] scheme. These features enable the model to capture local scale flows over complex terrain at high resolution. The model uses the second-order finite-difference method of Arakawa [1966] for treating the momentum equations in a way that conserves both momentum and kinetic energy. Thermodynamic and scalar conservation equations use the positive-definite advection scheme of Smolarkiewicz [1983, 1984] along with a non-oscillatory correction to suppress unphysical over- and under-shooting. In the model interior, a sixth-order spatial filter is used (with fourth-order and second-order filters active at the boundaries), which strongly damps motions with spatial scales less than approximately 4Ax, 4Ay. The M P I environment has been incorporated into WFIS, allowing multi-processor computation on distributed memory computers. In the M P I framework, different layers may have different resolutions whereas each sub-domain of any layer will have identical grid resolution in the three spatial dimensions, because it is the sum of the respective subdomains that constitute a layer. In this study, the numerical experiments were performed on the U B C Geo-Disaster Center I B M cluster, which has 256 processors, nicknamed "Monster". We use 16 processors for our numerical simulations. The M P I environment greatly improves the computational efficiency of the WFIS model. The model has been widely used to investigate and simulate local-scale atmospheric mechanisms, such as the simulations of deep convective cloud [Clark, 1982; Grabowski et al., 1991; Lane et al., 2003]; severe down-slope wind storm [Peltier and Clark, 1979; 1983; Clark and Farley, 1984], clear air turbulence  [Clark, 2000], weather modification  [Redelsperger and Clark, 1990; Hoinka and Clark, 1991], and wildfires [Clark et al, 2003]. These studies have demonstrated the ability of this model to simulate local-scale dynamics with reasonable accuracy. To simulate canopy flow at the West Coast Flux station, an important issue is the treatment of tree drag induced by the canopy elements, and the thermal energy budget  within the canopy and at the surface. The inclusion of additional physical parameterizations to describe these effects is an important part of this work. The remainder of this chapter focuses on the treatments of canopy dynamics and thermodynamics.  2.2 Treatment of canopy drag The canopy disturbs the airflow due to tree-drag force generated by the friction between the leaf and the airflow. The crucial difference between the flow inside the canopy and that above it is that considerable momentum is absorbed over the depth of the canopy in the form of drag due to the canopy elements. Flow above a canopy (z>1.5/z ) has been 0  observed to follow the log law or the modified log law depending on the atmospheric stability. Just above and in the upper canopy, however, the observed profiles deviate systematically from the theoretical (modified log law) profile [Arya, 1988]. Mean wind speeds within a forest canopy are well known to be considerably smaller than those over trees due mainly to the drag induced by leaves, stems and branches [Oliver, 1971; Zeng and Takahasshi, 2000]. Amiro [1990] stated that the velocities decrease more quickly within a denser canopy, and velocity distributions are skewed and kurtotic within the canopy. A l l these unique features in the canopy wind field are considered to result from the tree-drag force. In this section we first introduce a number of new variables, which are added to the model to describe the structure of canopy elements; and then describe the parameterization schemes used to treat the dynamical forces induced by the canopy elements in the model. The basic dynamical equations in the model related to this research, such as the momentum and mass conservation equations in terrain following coordinates, are described in Appendix A .  2.2.1 Basic variables describing canopy To simulate the flow above and within the canopy at the West Coast Flux station, a detailed description of the physical environment is needed. The description should span the range of scale from topographical features such as hills and coastlines down to small-scale horizontal and vertical distribution of canopy elements, because all these scales can affect the flow field at a given point such as flux tower site. The topographical features will be  described in detail in Chapter 3 along with a description of the experimental design. Here we focus on the structure of the canopy. The size and shape of canopy elements, as well as their distributions in space and time, determine the physical characteristics of canopy structure. However, including detailed canopy architecture information into numerical models is still in its infancy for two reasons: limited data on the structure of canopy elements and computational expense. A detailed description of non-homogeneous or non-uniform plant canopies in complex terrain would necessarily involve many parameters. A commonly accepted simplification is made by assuming that the plant stand is horizontally uniform and its average characteristics may vary only with height above the ground. In the model, the canopy structures are represented by the following new variables: •  The average height of a canopy, h, (33 m for the mature stand at the West Coast  Flux station). •  The cumulative leaf area index (LAI, units: m leaf/ m ground), defined as the one2  2  sided area of the leaves (needles) within a vertical cylinder of unit cross-section and canopy height h. It is the simplest useful measure of the canopy area density. LAI equals to 7 for the stands in our model. •  The foliage area density function, (a(z), units: m m" ), a local characteristic of the 2  3  canopy architecture, which represents the one-sided leaf area per unit volume of the canopy at height z. The relationship between LAI and a(z) is [Ross, 1975]: LAI=\a(z)dz  (2-1)  0  •  The local leaf area index, L(z), units: m leaf/ m ground, a more convenient local 2  2  characteristic, represents the leaf area per unit horizontal area of the canopy above the height z. h L{z) = ja(z)dz  (2-2)  Note that L{0) = LAI. •  Dry mass of the canopy needles/leaves and twigs per unit ground area, Com, (kg/m ),  a variable which characterizes canopy material.  2  •  Moisture content of the canopy needles/leaves and twigs per unit ground area, CMM,  (kg-m~ ), another variable to describe a canopy property. Combined with the dry mass, the 2  canopy density p (total leaf mass/ground area) can be determined. c  •  The fraction of ground area covered with trees, r|, an adjustable parameter, enables  the model to treat various cases such as a dense forest or a sparse forest or bare soil. r\ = 0 represents the bare soil and r\ — 1 represents a dense canopy. •  Canopy drag coefficient for an individual canopy element, C„d, is a dimensionless  parameter in the tree-drag formula.  2.2.2 Treatment of canopy drag Mathematically, Plate [1971] suggested the drag force, F, acting on an individual canopy element of effective area A per unit volume, and exposed to wind velocity, V, is (2-3a)  F = ± C A\V\V Pa  where p  a  nd  is the density of air, and A is the effective aerodynamic surface area of the  vegetation per unit volume. In this study, the drag F imposed by the canopy at height z within a canopy follows a volume drag formulation similar to Wilson et al. [1998] as F {z) = -C a{z)\V\v {z) i  nd  (2-3b)  i  where IV I is the absolute wind speed, calculated as ( Z v ; V , ) , / denotes the x, y, z 1/2  directions. The negative sign assures the direction of drag force is opposite to the wind direction. This parameterization scheme has been widely used in studies on wind flow within a vegetation canopy [Meyers and Paw U , 1986; Wang and Takle, 1996]. The tree drag appearing in the current model is: F {z) = -i C a{z)\V\v {z) i  1  nd  i  (2-3C)  The tree drag (2-3c) is implicitly treated in the momentum equations [(A-l), (A-2), (A-3) in Appendix A ] , acting as a momentum sink. For the Mature Stand site, the value of LAI we used was 7 during its full leaf seasons, according to the data from U B C Biometeorology/soil physics group. The vertical distribution of leaf area density a(z) in (2-3c) is parameterized as:  a(z)=(7/5)(l-(z-23)/10)  for 23<z<33  (2-4)  The cumulative leaf area index (the vertical integration of a{z)) was 7 from the canopy base (23 m) to the canopy top (33 m), i.e. LAI = \a(z)dz = 7. The vertical leaf area distribution 23 is shown in Figure 2.1. 33r 30 25 23 20  Its-4—'  -£Z O) 10 5 Q\  -1  ^  -0.5  .  0  .  0.5  1  .  1.4  Leaf area density (m /m ) 2  3  Figure 2.1: Vertical distribution of leaf area density (m" ). Vertical integration produces a leaf area index (LAI) of7m /m . 1  2  2  A continuing problem in the specification of the drag coefficient for an element embedded within a foliage canopy is that the element usually lies within the wind shadow of upwind elements, and its effectiveness as a drag elements may be substantially reduced, namely, the shelter factor. Thus, estimates of element drag coefficients based on calculations of the momentum balance within canopies yield values that are considerably smaller than those deduced from exposure of individual elements in wind tunnels [Thom, 1971]. The effective C„d of a single leaf measured in the wind tunnel also changes with leaf orientation and turbulent scales and intensity around the leaf [Raupach and Thom, 1981]. Currently, practical models used to simulate canopy flow are all based on appropriate assumptions for C„d. For example, Wilson et al. [1998] suggest using C„d = 1, while Shaw et al. [1988] used an isotropic drag coefficient of 0.15. However, Meyers and Paw U [1986]  obtained a constant value of C„d for their canopy by matching the model prediction of mean wind against observed canopy winds while Wilson [1988] obtained C„</ from the observed shear stress divergence. Thus, C d in the first case is used as an adjustable model parameter, and in the second as n  a data input. This points to one of the major weaknesses in our current understanding of canopy flows: we cannot predict the drag on the foliage from knowledge of the canopy geometry and the behavior of the plant elements in isolation. Brunet et al. [1994] have discussed the problems of canopy drag in detail. To bypass this difficulty, one can look at the sensitivity of varying C„d over a range of values to assess the effects on the dynamics for different flow regimes. Therefore, after we treated tree drag in the momentum equation, we performed the numerical experiments of sensitivity to drag-coefficient (C„d), which are described in Appendix C.  2.2 Refinement of model physics The description of the first law of thermodynamics and long-wave radiative-cooling schemes in the model are included in Appendix B . The model code contains several options to calculate long-wave radiation for no canopy. The radiation-cooling scheme used in this study is the combined R-C scheme [Redelsperger and Clark, 1990] (details in Appendix B), which accounts for the effects of carbon dioxide, water vapor and liquid water, and also is computationally efficient. The presence of a forest canopy significantly affects the surface energy budget. It shades the ground and traps heat, modifying air temperatures within the canopy and the ground. The R-C radiation scheme must be modified to include the effects of the canopy on radiation at the surface. In this section we describe how the radiation effects of canopy elements are incorporated into the surface and canopy energy budgets.  2.3.1 Energy budget within a canopy Radiative transfer in a plant canopy is a very complicated problem for which no satisfactory general solution has been found. The reason is that there are significant amounts of internal radiative absorption, reflection, transmission, and emission within the canopy elements, which vary with canopy architecture [Oke, 1987].  The parameterization scheme we choose to treat the canopy radiation effects is Yamada's [1982]. The net radiation R i, at treetop (at height h) is given by: N  R  Nh  =(l-a )S  +  l  £ c  (R  L h l  -R )  (2-5)  m  where a is the canopy albedo (= 0.1), and here S represents the direct solar radiation flux t  density at the treetop, the diffuse short-wave radiation is neglected under fair-weather for simplicity. WFIS model also has other options (like C C M 2 radiation code) for calculation of solar radiation under cloudy situation. e is the emissivity of the trees (= 0.98), c  R and Lhi  R f denote the incoming and outgoing long-wave radiations, respectively, which are Lh  calculated by using the R - C radiation scheme (Appendix B). The direct solar radiation S at the treetop in (2-5) is computed from: S= S  where S  m  cos Z  m  (2-6)  is the direct beam solar radiation at the treetop, and Z is the zenith angle of the  sun and is determined from the following formula: cosZ = sin^sin£ + cos^cosJcos//i?  (2-7)  where <f> is the latitude, S is the declination of the sun, and BR is the solar hour angle, positive clockwise from apparent noon, which depends on the model initial time. The direct beam solar radiation at treetop in (2-6) is specified as: S =S G m  (2-8)  Q  where So is the incoming radiation flux at the top of the atmosphere (= 1368 W m" ). The 2  fraction (G) of solar radiation reaching the surface is much less than that at the top of the atmosphere due to many factors, including molecular scattering and absorption by permanent gases such as oxygen, ozone, and carbon dioxide. This effect is parameterized by Atwater and Brown (1974) to include the effect of the forward Reyleigh scattering. The fraction is: G = 0.485 + 0.515 1.041 -0.16( where P  R S  0.000949P,., +0.051 cosZ  1 / ? 1 2  )'  (2-9)  is the surface pressure in millibars.  Combining (2-9), (2-8), (2-7), (2-6) and the R - C long-wave radiation scheme, the net radiation R h at the treetop can be calculated. N  The net radiation RN within a canopy is following Yamada [1982], based on Uchijima P  [1961] as quoted in Ross [1975]: R {z) Np  = ^[exp{-^(z)}-7(l-|)exp{-Ai(0)}]  (2-10)  where r\ is a fraction of ground area covered with trees, k is the extinction coefficient (= 0.6 in this study as in Yamada (1982)), RNH is calculated using (2-5), and h is the canopy height. L(z) is the local cumulative leaf area index, which represents the leaf area per unit horizontal area of the canopy above the height z, L(0) = LAI = 7. Equation (2-10) states that the transmission of incoming net radiation into a plant canopy shows an approximately exponential decay with the depth of penetration, as attenuated by the depth-dependent local leaf area index. The energy budget within a plant canopy is: R =H  + LE + S,+P  Np  (2-11)  s  where RN is the net radiation within the canopy, H is the sensible heat flux, LE is the latent P  heat flux, S is the heat storage in the canopy, P is the energy consumed via photosynthesis. T  S  Following Yamada [1982] and Song et al. [1997], as a first approximation, the storage term S is neglected by assuming that all the energy stored in the stems and branches T  becomes available immediately to heat adjacent air. The energy used in photosynthesis is less than a few percent of the net radiation [Oke, 1987; Twine, 2000] and is neglected here for simplicity. Under the conditions assumed, the internal energy equation within a canopy is written as: ^ ^  =  ( 1 ^ ) ^ Pa p C  8 2  L  1  +  :  Pa p+Pc c C  C  ( 1+  I .^aL r  B  ( 2  .  1 2 )  <%  where the long wave radiation RN without canopy is calculated by using R - C long-wave radiation scheme (details in Appendix B), and R  NP  within the canopy is calculated by (2-10).  The Bowen Ratio B (= HILE) within a canopy is assumed to be constant with height, as is supported by the observations of Bradley et al. [1980]. In our simulations, the Bowen Ratio is set to 1.5 within the canopy [Yamada, 1982], and over water is set to. 0.1. p and p a  c  represent air density and canopy density, respectively; C and C is the specific heat of air p  c  and canopy, respectively. C / C is set to 0.65 currently in the model, r\ is a fraction of an c  p  area covered with trees, r\ = 1 represents a dense forest, and r\ = 0 for bare soil. The local canopy density p (z) is related to the leaf area index and to the canopy dry and c  moisture masses, and is parameterized as p (z) = a(z)*(C c  +C )  Dm  (2-13)  Mm  where a(z) is volume leaf area of the canopy. C , dry mass of the canopy needles/leaves Dm  and twigs, is set to 2.0 kg/m . CM™, moisture content of the canopy needles/leaves and twigs, 2  is set to 2.5 times of dry mass, i.e. 5.0 kg/m in our simulations. 2  2.3.2 Soil surface heat budget Exchanges of energy between the atmosphere and the earth's surface  fundamentally  influence the dynamics of the atmosphere. Unfortunately, we cannot use the same parameterizations  as within the canopy because the heat conduction is much more  significant between a solid surface (soil) and its neighboring air. The surface energy balance in the current model is treated as: R =H +LE NG  C  G  +  G,  (2-14)  RNG is the net radiation to the ground, which comprises net short wave and net long-wave radiation, H  G  and LEQ are the sensible and latent heat fluxes to (or from) the atmosphere,  and G represents the heat flux into or out of the sub-medium (soil or water). Equation (2-14) S  describes how the net radiation at the surface must be balanced by a combination of the sensible and latent heat fluxes to the air and the heat flux to the subsurface medium (soil). Considering the canopy effect on the net radiation budget at the surface, RNG in (2-14) is given by: RNG = vR  Nh  exp[-*£(0)] + (1 - / ? ) [ ( ! - a ) S + s (R G  G  LGL  -R  L C N  )]  (2-15)  where « c is the surface albedo (set to 0.25 for water, 0.3 for bare ground), S is the direct solar radiation (2-6); R  LCI  represents the incoming long-wave radiation to the ground, R ^ LC  is the outgoing long-wave radiation from the ground; 8 q denotes the infrared radiation emissivity of the ground (set to 0.98 in the model). L(0) = LAI = 7. In (2-14), the sensible heat flux HQ and latent heat flux LEQ between the ground and the adjacent air can be expressed as by Devonec and Barros [2002]:  H =C K \U \(T -T )  (2-16)  LE =L K \U,\{q,-q ,)  (2-17)  G  c  pPa  H  vPa  x  x  0  y)  sa  The subscript 1 indicates a reference height in the boundary layer at which the horizontal wind U\, the temperature T\, and the specific humidity q\ are known. In this study, we use the values of the first layer above the ground. Ku and A"w are the aerodynamic drag transfer coefficient for heat and moisture (dimensionless), respectively. Currently we assume their values are the same. g t is the saturated specific humidity at the ground surface, Ly sa  represents the latent heat for water vaporization or condensation, and Cp denotes the specific heat of air at constant pressure (units: J-kg"'-K~ ). To is the ground surface l  temperature. The rate of heat transfer or heat flux in a given direction is found to be proportional to the temperature gradient in that direction, i.e., the heat flux in the z direction, Gs=-C k ^s  dz  s  (2-18) s  where z is the depth of the soil and positive downward (in meters), k is the soil thermal s  s  diffusivity ( m V ) , and C is the vo^ volumetric heat capacity of the soil (units: J-m" -K ). 1  3  _1  s  Following Yamada [1982], we adopt: k = (0.\5 + 4z )*10~  6  s  s  (2-19)  The soil temperature T is obtained by s  ^=!-(*. I) 1  -  (2 20)  The boundary conditions for (2-20) are the heat energy balance at the soil surface and the specification of a very weak constant heat flux of soil at 77 cm depth in the soil. (2-20) is numerically used in 4 soil layers: 0 - 5 cm (5 cm thick), 5 -17 cm (12 cm thick), 17 - 41 cm (24 cm thick) and 41 - 77 cm (36 cm thick).  2.4 Summary The presence of a canopy leads to a number of numerical issues in the simulation of boundary-layer flow in a deep forest. Canopy elements absorb momentum, trap and consume solar radiation and shade the ground. These cause significant changes to the momentum and energy budgets within the canopy, and eventually affect the profiles of the  flow and scalar fields within and just above the canopy. However, the original code of WFIS didn't include the parameterization schemes describing canopy structure and its dynamic and thermodynamic effects on the boundary layer flow. Thus, part of the work was to add these effects to WFIS. After considering a number of parameterization schemes, we deal with tree drag in the momentum equations following Wilson [1998], treat the energy budget within the canopy following Yamada [1982], and deal with surface heat budget in the manner suggested by Devonec and Barros [2002]. To balance the surface energy budget and get ready for the treatment of a CO2 source term in the CO2 conservation equation, we incorporate a fourlayer soil model using a heat transfer equation to simulate the soil temperatures at different depths. In the course of the refinement of the model, we also performed a series of numerical sensitivity experiments to test the model code and to examine the effects of some isolated forcing processes on canopy flow. The preliminary results of these numerical tests are included in the Appendix C.  Chapter 3 Simulations of canopyflowover a hill 3.1 Introduction Wind flow is an important factor for transporting scalars (heat, water vapor, carbon dioxide, pollutants, etc) within a forested canopy. Understanding the mean and turbulent wind regimes at the West Coast Flux station is a prerequisite for interpreting the flux measurements at the tower site. Most investigations of flow inside and above plant canopies have been experimental, including in-situ field measurements or wind-tunnel experiments [Seginer et al., 1976, Finnigan, 1979; Shaw et al., 1988; Raupach et al., 1991; Ruck and Adams, 1991]. The wind profiles above plant canopy (z > 1.5h) are observed to follow a log law or modified log law, depending on the atmospheric stability [Stall, 1988; Kaimal and Finnigan, 1994]. However, just above and within the canopy, the observed wind profiles deviate systematically from the theoretical (modified log law) profile [Arya, 1988]. The flow inside a canopy is rich in features arising from the interactions between the canopy and the atmosphere, exhibiting high variability both spatially and temporally, as shown by almost two decades of published research [e.g. Raupach et al., 1986; Gao et a l , 1989; Amiro, 1990; Wang and Takle, 1996; Wilson et al., 1998; Novak et a l , 2000; Sakai et al., 2001]. To describe mean sub-canopy wind profiles, it has traditionally been assumed that the forest is horizontally homogeneous and that in the denser canopies (LAI > 1) all the momentum is absorbed in the upper part of the canopy, so the shear stress transmitted to the ground is nearly zero. The upper part of the wind profile within the canopy can be approximated by: V(z)=V e- ^ r  h  (3-1)  where V(z) represents the wind speed at height z, Vi, is the wind speed at the top of the canopy, h is the height o f the canopy, and T is an extinction coefficient, depending on the canopy leaf area density and the mixing length in the canopy. Inoue [1963] and Cionco [1965] also derived the wind profile under the assumption of a constant mixing length.  However, Kaimal and Finnigan [1994] summarized wind profiles normalized by the wind velocity at canopy height for a variety of real and wind-tunnel canopies, and found that the normalized wind profiles do not simply collapse to a common curve, especially for a dense canopy. They attribute this to significant vertical variation in foliage density in some real canopies. Accurate simulation and prediction of canopy flow has been handicapped due to the complexity in the array of vegetation elements (leaves, branches, etc.) and the complex processes of air momentum transport within a vegetation canopy. In the case of the forest canopy located on complex terrain, the situation is even more complicated. There is very little in the literature on the impacts of tall trees on air circulations over sloping terrain. Some studies discussed the effects of idealized hills on the flow fields at different positions along the slope, but were restricted to bare hills [Taylor et al., 1987; Raupach and Finnigan, 1997]. Others tried to simulate turbulent flow within the canopy with various closure models, but they focused only on a homogenous canopy over flat terrain [e.g. Wilson and Shaw, 1977; L i et al., 1985; Wilson, 1988; Miller et al., 1991; Ayotte et al,, 1999; Finnigan, 2000]. The studies of flow over a ridge covered with a deep plant canopy are limited to wind tunnel experiments [Ruck and Adams, 1991; Finnigan and Brunt, 1995], but there are problems in reproducing realistic atmospheric stability conditions. Only recently have some publications addressed canopy flows over forested hills. Based on Hunt et al. [1988a; b], Finnigan and Belcher [2002] developed a linear analytical model for a neutrally stratified A B L flow over a hill covered with a homogeneous canopy. Their results indicate that the maximum winds are near the crest of the hill. In the lee of the hill, an adverse pressure gradient leads to flow separation, with a region of reversed flow when the canopy is deep enough. Brown [2001] and Allen and Brown [2002] conducted largeeddy simulation of neutral turbulent flow over two-dimensional ridges. Issues related to resolving the flow in the dynamically important inner region led them to represent the roughness elements on the surface of the hill using a canopy. They found that the canopy reduced the speed-up of the flow over the hill, and increased both the drag on the hill and the tendency for the flow to separate. Therefore, although the airflow inside and above plant canopies has been studied intensely during the past decades, our understanding of atmospheric boundary layer ( A B L )  flow over plant canopies is restricted to flow over homogenous canopies on level terrain. There is a fundamental need to better understand canopy flow over complex terrain. Further, parameterization of topographic roughness has been shown not only to improve the skill of large-scale numerical weather prediction models [Belcher and Hunt, 1998], but also to make advances in related applied sciences, such as air pollution, hydrology, and forestry. This study employs a high-resolution large-eddy simulation model, the Weather Fire Integrated System (WFIS), to investigate canopy flow over complex terrain. With finer grid resolution, the averaging volume is sufficiently small so that the larger, energy-containing scales of motions are explicitly resolved, with the beneficial result that the contribution to the dynamics by the sub-grid-scale eddy-mixing parameterization is small. A l l numerical experiments in this research are performed in two-dimensions (2D). This framework is useful for understanding the complex effects of tall trees on air circulations in the atmospheric boundary layer over sloping terrain, and it provides theoretical and physical guidance for more computationally intensive three-dimensional calculations in the future. Also, in two dimensions, higher resolution can be attained in comparison to threedimensional (3D) calculations. In the course of the refinement of the model, we also conducted a series of numerical sensitivity experiments to evaluate and to diagnose the performance of the parameterization schemes described in Chapter 2. These preliminary numerical experiments examined the response of canopy flow to some isolated processes such as the obstacle effects of topography on the ambient flow, the impacts of an isolated canopy on the lee flow, and the response of ambient winds to tree drag. These sensitivity experiments and their results are described in the Appendix C. The results of these sensitivity experiments suggested that an isolated canopy on the lee side of a hill increases the tendency for the lee flow to separate and reverse i f the tree-drag is sufficiently large, i.e., for a sufficiently dense canopy. The tests of sensitivity to upstream flow and topography indicate that higher hills result in stronger flow impinging on the canopy, and lead to earlier flow separation and stronger flow reversals at the downstream edge of the canopy. The interaction between the mountain and the canopy makes the properties of the flow more interesting and also more difficult to study. To interpret  measurements made near hills and windbreaks, one has to consider all these factors and distinguish which factor is dominant before reasonable explanation can be made. This chapter describes numerical experiments designed to examine the influence of a dense canopy on the surface energy budget and on the A B L flows over sloping terrain in the context of a diurnal cycle. We will gain insights on the local-scale controlling mechanisms that are superimposed on the photosynthetic and respiratory processes of the ecosystem, and which might account for the anomalies in the E C flux measurements at the West Coast Flux station.  3.2 Experimental design The flux measurements are affected by myriad transport mechanisms, which vary on numerous spatial and temporal scales. The dominant scale is the diurnal cycle with generally statistically unstable conditions during the day and stable stratification at night. The presence of a deep canopy could have significant impacts on the surface energy budget because the canopy elements absorb momentum, trap and consume solar radiation, and shade the ground. Two series of experiments were carried out to simulate the mesoscale A B L flows at the West Coast Flux station under a diurnal context. Series one focused on the effects of a tall canopy on the surface heat budget, with the dynamics turned off to allow a longer time step and longer model-run duration. Series two was a full two-dimensional application of the WFIS code with both dynamics and physics included, to investigate the diurnal and topographic response of canopy flow. Diurnal variation of all variables was computed by taking the canopy radiation balance into consideration. These experiments outlined in Table 3.1 employed a two-way interactive nested domain. Non-standard Cartesian notation was used in the WFIS runs, resulting in y representing the eastward direction rather than x. The domain configurations are shown in Fig. 3.1, which mimic a west-to-east cross-section of Vancouver Island through the Mature Stand site at the West Coast Flux station. The horizontal size of the outer domain was 86.4 km, within which the land surface extended from y = 11 km in the west to y = 64 km, and which contained an idealized mountain range of height 1300 m, representing Vancouver Island. A canopy of trees 33 m high was also specified. The leaf area density, a(z), was specified  using (2-4). On the either side of the land was water: from y = 0 km to 11 km was the eastern edge of the Pacific, and from 64 km to 86.4 km was Georgia Strait. The finer, nested domain was located on the east side of the mountain, simulating the location of the flux site. A 1 km deep boundary layer and calm wind was assumed, and the atmosphere was taken as weakly stable (Brunt-Vaisala frequency, N = 0.002 s" ) from z = 0 to 1 km. The 1  free atmosphere was assumed to be N = 0.01 s" from z = 1 to 13.5 km, and 7Y = 0.02 s" 1  1  above 13.5 km. The Coriolis force was turned off. Because the anomalies in the actual flux measurements show up during fair weather, our numerical experiments assumed clear sky. At z = 0 an initial surface temperature of 288 K for the land and 287.5 K for the water was specified.  37.8  46.4  55.1 Y  63.7  72.4  81  ( KM )  Figure 3.1: Schematic diagram of topography and canopy specifications for the outer domain (upper) and the nested domain (lower), where the shaded area represents the canopy. The nested domain is located on the eastern side of the mountain range. Table 3.1 outlines the general characteristics of the diurnal-cycle experiments. Ay and Az represent the grid size in the horizontal and vertical, respectively. A? denotes the time step, r| represents the fraction of ground area covered with trees, and is an adjustable parameter enabling the model to treat various cases such as a dense forest or a sparse forest or bare soil. In these experiments, r\ = 0 represents the bare soil, r\ = 0.8 denotes a sparse canopy and T) = 1 represents a dense canopy. A leaf area index of 7 is specified. Duration represents the model simulation time.  Although there are some arguments on the effective drag coefficient, the preceding sensitivity experiments to drag coefficient (see Appendix C) suggested that a drag coefficient between 0.1 and 1.0 made only minor difference to the disturbed flows. With this in mind, we chose C„d = 0.15 for the following experiments outlined in Table 3.1. This value was also suggested by Amiro [1990], and is considered typical for mixed deciduous forests. The dynamics was turned off for those experiments labeled as N O F experiments, which were run for 100 hours duration starting at 0500 L S T . These simulations assumed solar radiation valid for late July (July 30th). The purpose of the N O F experiments was to test the surface energy budget code, which included the effects of the canopy on the surface radiation budget, and included the effects of soil conduction. Experiments labeled as INF are the same as N O F , except the dynamics were active, and they were run for 42 hours duration. Table 3.1: outlines the characteristics of diurnal-cycle experiments. Ay and Az represent the grid size of horizontal and vertical, respectively. At denotes the time step. Slashes separate values for the two nested domains. Exp. Ay(m) Dynamics H i l l height Duration At (s) Az j 11 (m) (m) (hours) (tree-cover) NOF0 150/75 12/6 20/10 0 off 1300 100 NOF1 12/6 100 150/75 20/10 0.8 off 1300 NOF2 150/75 12/6 20/10 1 off 1300 100 INFO 150/75 12/6 0 on 1300 42 1.0/0.5 INF1 42 150/75 12/6 3.0/1.5 0.8 on 1300 INF2 1 on 42 150/75 12/6 3.0/1.5 1300 m  n  3.3 Vertical grid design The model employs vertically stretched, terrain following coordinates with grid refinement. The broad range of temporal and spatial scales involved in simulating canopy flow in the lee of a hill appears to be a very challenging numerical problem. When we try to use higher resolution to evaluate the dynamic response of the mountain-canopy to some specific mountain flow regimes, strong shears induced by mountain drag and canopy drag lead to high frequency oscillations (see Appendix C), resulting from the second-order numerical treatment of momentum advection. One solution to this problem is to include a nonoscillatory advection correction for the momentum advection terms, but a non-oscillatory  correction is already included in the model for the advection of the thermodynamic variables through the use of the Smolarkiewicz [1984] and Grabowski and Clark [1990] advection scheme. Another approach to solve the problem is to redesign the vertical stretched grids. Clark and Hall [1996] stated that an important consideration in the design of vertically stretched grids, particularly with vertical grid refinement, is avoiding model artifacts  resulting from noisy vertical derivatives of the Jacobian of the vertical  transformation. The vertical layers interact with each other at the quadratic or cubic level, resulting in irregularities. When the initial grid is defined using a high-order hyperbolic tangent profile, it is important that the Jacobian satisfies the same conservation principles used in the model interpolation formula, which means that only the outermost domain can be independently designed.  20 .00  13 .33 I 6 .67  0 . 00 0.0  125.0 250.0  375.0  0 . 00 0.0 500. Z  10.0  20.0  30. e  42 <y) az <M) Figure 3.2: Vertical grid structure for the two nests, (a) Vertical grid size Az versus z for nest one, (b) same as (a) but for nest two.  The sensitivity experiments in Appendix C suggested that the vertical velocities associated with the canopy drag were strongest at an elevated position. Thus, the vertical grid increment should be designed so that the highest amplitude values of vertical velocity occur in a region of relative coarse vertical grid increment. Figure 3.2 shows the vertical grid structures employed.in the diurnal cycle experiments, which are derived by first specifying a smooth vertical grid for nest one. The inner nests are obtained by projecting the outer nest into the inner nests. The outer nest is first specified with 4 layers of constant Az = 12 m followed by 9 layers having linear increase to Az =150 m, then followed by 75 layers of Az =150 m followed by another 6 layers with linear increase to Az = 450 m, and finally followed by 16 layers of Az = 450 m as shown in Fig.  3.2 (a), which extends to over 21 km. The outer nest has 110 vertical layers in total. The inner nest starts from Az = 6 m and gradually increases to Az = 75 m, there are 30 vertical layers for the nested domain, which extends to 1.2 km high. Since the outer domain starts with Az = 12 m at the surface and quickly increases with height, the large amplitude of vertical velocities have a much less restrictive effect on the time step. This vertical grid design results in a broad range in vertical resolution suitable for the treatment of canopy flow.  3.4 Results and discussion 3.4.1 Surface temperature variations Figure 3.3 shows a comparison of the variations of the soil surface layer (2.5 cm thick) temperature and the air temperature within the trunk space (6 m above ground) for a fourday diurnal cycle simulation with dynamics off for experiments N O F [0-2]. Figure 3.4 is similar to Fig. 3.3 except it is a two-day simulation with dynamics active for experiments INF [0-2]. The temperature analysis is based on a single horizontal point (grid 405 located at y = 60.67 km). This calculation is representative of the conditions under a mature 33 m high canopy with LAI = 7 at the West Coast Flux station. For bare ground (NOF0, r\ = 0), the range of soil and air temperature is around 20 K and 16 K . For the denser canopy cases, the range is about 12 K for r| = 0.8 and about 9 K for r| = 1 (Fig. 3.3). A dense canopy reduces the diurnal range of surface temperature because a large part of the incoming solar radiation is absorbed by the plant surface, reducing the amount of solar radiation reaching the ground and in turn decreasing heating to the neighboring air. Hence, the surface temperature during the day under a plant canopy is uniformly lower than over a bare soil surface. At night the outgoing long-wave radiation is partly intercepted and absorbed by vegetation, which also radiates energy back to the surface. This slows down the radiative cooling of the soil. The combined effects of the two processes result in smaller diurnal range of temperatures on vegetative surfaces. The increasing tendency of temperature in Fig. 3.3 is because advection is turned off, allowing solar radiation to exceed the long-wave radiation loss, resulting in heat accumulation within the canopy.  320.0  320.0  = 0.8  313.3  313.3  A  300.0  "  293.3 286.7  J\ A /  55.0 TIME  IHRI  0 105.0  i  ,  ,  1  1 -I  i  30.0  ti=i.o  A  293.3 \  -  i—i—i—i—  55.0 TIME  !  ? 300.0  i  1 r- ,  '-  '  N  280.0 30.0  :  ,  306.7  306.7 J  ,  IHRI  if\ '• '  286.7 280.0 30.0  55.0 TIME  B0.0 105.  <HR)  Figure 3.3: Temperature of soil surface layer (2.5 cm, thick line) and surface air temperature (3 m, thin line) for NOF0 (r| = 0), NOF1 (r| = 0.8), and NOF2 (r) = 1). Dashed lines mark midnight. With dynamics active, the diurnal range of temperature about 12 K (r\ = 0), 6 K (r\ = 0.8) and 4 K (n = 1.0) is much smaller than that of N O F [0-2] with dynamics off (Fig. 3.4). There are two processes causing the difference. One is advection that transports heat away from the canopy space. Another is the enhanced turbulence by the canopy near the surface, providing more effective exchanges of sensible and latent heat between the soil and overlying air. The combined effect of the two processes significantly reduces the diurnal range of surface temperature for the case with dynamics on. The maximum temperature for a dense canopy occurs around one hour later than that for bare-soil case, and the minimum temperature is reached in the early morning when the net solar radiation balances the net long-wave radiation. The increase of temperature from one day to the next is much smaller, because of the effects of advection and turbulent mixing. Comparing the air temperature to that of soil surface layer in Fig. 3.3 and Fig. 3.4, it is evident that for bare soil (r\ = 0) the temperature of soil surface layer is cooler than air temperature at night, and is much warmer during daytime. This is because the soil is a better absorber of solar radiation and also a better emitter of long-wave radiation compared to the air. For the denser canopy (r\ = 0.8 and 1.0), the situation is much different. The surface air temperature is a little warmer than the surface layer of the soil during daytime, while at night they are almost the same. The difference is due to the radiation effects of a dense canopy on the surface energy budget. During daytime, most of incoming solar radiation is absorbed by the plant canopy, which in turn heats the surrounding air and causes the air just above the canopy to be the warmest. With strong turbulent mixing induced by the canopy elements, the air in the trunk space becomes warmer than the ground.  At night, the cold canopy top creates an unstable stratification below it, which leads to the vigorous eddy mixing within the canopy, and results in little temperature  difference  between the surface air and the surface layer of the soil. 310.0  n= 0.8  305.0 300.0 295.0 290.0 285.0 280.0 15.5  26.0 TIME (HRI  36.5  17.  15.5  26.0 TIME (HRI  36.5  47  Figure 3.4: Similar to Fig. 3.3 but for experiments INF [0-2] with dynamics on.  3.4.2. Katabatic and anabatic flow In addition to perturbing the background flows as we showed in the sensitivity experiments (see Appendix C), hills can also induce their own thermally driven topographic flows; namely, gravity flows. The West Coast Flux tower site is located on a 5-10 degree slope oriented east/northeast on Vancouver Island, where katabatic and anabatic winds induced by differential heating and cooling are expected to be strong during fair weather. Figure 3.5 shows the simulated horizontal winds at four times spaced 6 hours apart for flow over Vancouver Island for experiment INF2. The daytime convective anabatic winds and nighttime stable katabatic flows are evident. A t 0600 L S T weak anabatic winds first appear on the eastern side of the hill because of the earlier solar heating on this side. A t 1200 L S T the strong convective anabatic flows occur on both sides of the hill, while at midnight (0000 LST) the katabatic winds dominate. But at 1800 LST, near sunset, the wind pattern is very complex and is in a transition phase as it switches from upslope to downslope flow. The wind pattern near sunset is shown in Fig. 3.6. At 1730 L S T anabatic flows still dominate on both slopes. A t 1745 LST the upslope flow on the eastern slope becomes much weaker, and the anabatic winds on the western slope are still very strong and are causing flow onto the eastern slope. B y 1830 LST, after the transition, the down-slope flows occupy the entire eastern slope, and are relatively weak.  O  10  20  30  40  Y  O  10  20  30  50  40  Y  O  IO  20  30  10  20  30  70  80  50  60  70  80  60  70  BO  (km)  40  Y  O  60  (km)  50  (km)  40  Y  50  60  70  80  (km)  Figure 3.5: Horizontal wind at the four times of 0600, 1200, 1800 and 0000 LST for the outer nest for INF2. Contour interval is 1 m/s with solid showing positive values from west to east and dashed showing negative.  Y  (km)  Figure 3.6: Similar to Fig. 3.5 except for the four times of 1730, 1745, 1800 and 1830 LST.  The dynamical explanation for this behavior is that it is a combination of the shadowing effects on the eastern slope (earlier sunset) and west-slope dominant upslope winds extending over the ridge because of the stronger heating rates on the western slope. The elimination of sensible heating on the eastern slope then allows a very quick transition to down-slope winds over the entire eastern slope in about one hour after sunset. The rather rapid shift from upslope to down-slope winds just after sunset may be an important factor in helping to explain some of the CO2 flux anomalies just after sunset at the West Coast Flux tower.  37.8  46.4  55.1 63.7 Y (KM)  72.4  81  37.8  46.4  55.1 63.7 Y (KM)  72.4  81  Figure 3.7: Drainage flow (horizontal wind velocity at night) for experiments INFO (left, n = 0) and INF2 (right, n = 1). Contour interval is 1 m/s with solid representing flow from west to east, dashed showing opposite. The topography is shown at the top of frame. 3.4.3 Drainage flow Nocturnal drainage flow is common in summer when the combination of clear skies and low synoptic wind-speed is frequent [Mahrt, 1986; Stull, 1988]. Most studies of drainage flow considered landscapes with sparse vegetation or bare ground [Manins and Sawford, 1979; Mahrt, 1982; 1986; Fitzjarrald, 1986]. In such cases, radiative cooling generates the lowest temperatures near the surface, and the drainage flow is typically stably stratified with a maximum velocity close to the slope surface. The principal forces generating the drainage flow are the (negative) buoyancy associated with cooling close to the surface, the drag at the ground, and the rate of turbulent transfer of momentum at the upper boundary of the drainage flow.  The structure of the drainage flow is more complex when slopes are covered with dense tall vegetation. The canopy top is the strongest source of radiative cooling to the night sky. As air cools near the canopy top it generates a stably stratified drainage flow above the canopy, similar to that on bare slopes but displaced upwards. The canopy provides frictional drag to this flow and generates turbulence that mixes air into the trunk space, while strong thermal stability associated with the canopy-top inversion inhibits turbulence persistence. The vertical gradient of air density associated with cold temperatures at the canopy top allows air to subside into the canopy. Down-slope gravity flow within the canopy is therefore likely, but would be much different from that over bare slopes. To give some insight on the nocturnal flows at the West Coast Flux station, the comparison of nocturnal drainage flows for INFO and INF2 is given in Fig. 3.7. For the bare soil, the nocturnal drainage flow is much stronger than that within the dense canopy, and extends to almost 200 m high with the maximum values of 6 m/s near the ground. When slope is covered by dense forest, the drainage flow is much weaker and shallower (1~2 m/s), and is isolated from the flow above the canopy by the temperature inversion at the canopy top. A n above-canopy drainage flow and a sub-canopy drainage flow appear. Recently, Komatsu et al. [2003] reported that this form of drainage flow was quite common over a tropical forest when wind speeds were low at night. Some studies indicate that flux measurements made over sloping terrain underestimate the nocturnal respiration on calm nights [Lavigne et al., 1997; Mahrt, 1998]. The sub-canopy drainage flow might be one of mechanisms that account for the transport of the soil-respired CO2 from land to the nearby water surface at the West Coast Flux station.  3.4.4 Thermals and Land/Sea breezes Figure 3.8 shows thermal structure and wind pattern at 1400 and 0200 local standard times (LST) for the nested domain for experiment INF2 (r\ = 1.0). A t midday (1400 LST) there are strong thermal updrafts associated with strong convection over the land, and sea breezes are evident. The sea breezes develop during daylight' hours, and reach a peak in midafternoon, and decay in the evening (also shown in Fig. 3.5). At night, light winds indicate a relatively stable boundary layer, where the land cools down much faster than the sea, and  weak land breezes form. The land/sea breezes might further amplify the up- and downslope flows, which transport CO2 away from the tower site.  -~  Y (km)  40  *5  SO  55  60  65  70  m/s  75  80  Y (km) m/s Figure 3.8: Thermal structure and land-sea breezes at 1400 and 0200 LST for nested domain for experiments INF2. Solid contours represent air potential temperature; contour interval is 1 K. The arrows denote winds. Shaded area represents the canopy.  40  45  50  55  60  Y (km)  65  70  75  80  m/s  Figure 3.9: Wind field at 1400 LST for the lower 150 m of Fig.3.8. Although anabatic flows and sea breezes dominate the flow pattern of the A B L during the daytime (Fig. 3.8), the wind vectors deep within the canopy suggest that the down-slope flow occurs frequently even at midday. To clearly see the flow within the canopy, Fig. 3.9 shows the wind fields at 1400 L S T for the lower 150 m, corresponding to the bottom portion of Fig. 3.8. From the scale ( 2 - 4 km) of the convergence area within the canopy, we infer that the intermittent downslope flows at midday in our simulations result from convective thermals turning over. The observations at the West Coast Flux station also indicate that above the canopy the wind pattern under fair weather displays upslope flow during the day and downslope at night. However, deep within the canopy the downslope flows dominate 75% of the time during 24 hours. To study such an effect, a further statistical analysis of the up- and down-slope flow over a 24-hour period inside the canopy is done next. Figure 3.10 shows time-line statistics of horizontal (a) and vertical (b) velocity at 6 m above the ground, horizontal (c) and vertical (d) velocity at 42 m above the ground at grid point 280 (located aty = 58.8 km in the nested domain), which has the same elevation 320 m as the Mature Stand tower site at the West Coast Flux station. During most of the daytime, horizontal velocities are negative (westerly wind) and vertical velocities are positive (upward), indicating upslope flows. Vertical velocity is well correlated with horizontal velocity.  15  20  29  Local time (hr)  Figure 3.10: Statistics of horizontal (a) and vertical (b) wind velocity at 6 m above ground, and horizontal (c) and vertical (d) at 42 m above the ground at grid 280 for experiments INF2 (n = 1.0). In (a) and (c) positive represents westerly wind (down-slope) and negative denotes easterly wind (upslope). In (b) and (d) positive represents upward (upslope) and negative denotes downward (down-slope).  There are mainly three periods during which down-slope flow occurs at daytime. They are around sunrise (0500 to 0700 LST), in the midday (1200-1400 LST) and at sunset (1700-1900 L S T ) . Thus, during daytime, down-slope flows usually appear at transition times and when strong convective thermals and turbulent mixing dominates. A t night the wind is light and dominated by down-slope flows. We output data from the model every 30 seconds, giving a total of 2880 data points for the 24-hour statistics in Fig 3.10. The statistics show that within the canopy down-slope flows occur 1725 times, while upslope flows appear 1155 times. Hence, during 24-hour period, down-slope flow dominated 60% of the time inside the canopy for this simulation.  3.4.5 Vertical profiles Vertical profiles of potential temperature show how the tall canopy modifies air temperatures within and above the canopy (Fig. 3.11). The daytime profiles (left panels in Fig. 3.11) indicate radiative warming of the canopy top resulting in a convectively unstable layer from the canopy top to 800 m high (Fig. 3.12). The lack of sunlight reaching the floor results in the colder air near the ground and a relative stable surface layer within the canopy even at midday. At night (right panels in Fig. 3.11), above the canopy top is a typical stable nocturnal boundary layer because the rate of radiation cooling at the canopy top is the biggest. Close to the surface the stability structure is reversed, the ground is a little warmer than the air temperature just above the ground because the canopy traps and also emits long-wave radiation back to the ground. The simulated horizontal and vertical wind profiles at grid 280 (same grid as Fig. 3.10, 3.11) for 1300 and 0100 L S T are shown in Fig. 3.12 and Fig. 3.13, respectively. The upper panels show the winds from the ground to 1000 m, and the lower panels show winds from the ground to 50 m. Above the canopy, the wind gradually increases, with maximum value occurring at 350 m altitude during daytime, and at 150 m at night, indicating a nocturnal jet. Inside the canopy, the minimum wind occurs in the upper canopy because of the large momentum absorption by the dense cluster of leaves. Farther down, as the leaf area density gradually decreases to zero, the tree drag disappears, allowing the buoyancy term to become more significant and dominate over the flux-divergence term. Consequently, the wind speed increases with decreasing height above the ground, until surface friction  reverses this trend close to the ground level. Therefore, a low-level wind maximum appears in the trunk space below the forest canopy (lower panels in Fig. 3.12, 3.13). Right above the ground, a weak down-slope flow appears both day and night.  potential temperature (K) Figure 3.11: Potential temperature profiles for daytime (left) and nighttime (right). The lower panels show air potential temperature from the ground to 50 m. In summary, during daytime, the surface temperature under a dense canopy is uniformly lower than that over bare soil, and the range of surface temperature for the dense canopy is much smaller. Above the canopy anabatic flows dominate during the day because of the effects of the sloping terrain. But near and within the canopy, the radiative effects of a tall dense canopy and strong thermal convection in the upper canopy transport warm air down into the lower part of the canopy, which induces intermittent down-slope flows even at daytime, but which mostly occurs at the transition times and when negative buoyancy  dominates. The statistics of simulated winds 6 m above ground indicate that 60% of time during 24-hours, down-slope flow occurs within the canopy. A t night, drainage flows dominate, but the drainage flow for dense canopy is quite different from that for bare ground. The dense canopy separates the drainage flow into sub-canopy (trunk space) and above-canopy parts, which is much weaker than that for the bare soil case. The heterogeneous surface (land and sea) results in land and sea breezes, which could further amplify the up- and down-slope flow. The minimum winds occur in the upper canopy because of the impacts of tree drag.  -6  -4  -2  0  2  0  2  4  Horizontal velocity (m/s) Figure 3.12: Similar to Fig. 3.11 except for horizontal wind.  6  Vertical velocity (m/s)  Figure 3.13: Similar to Fig. 3.12 except for vertical velocity.  3.5 Conclusions To understand the behavior of CO2 fluxes measured at the West Coast Flux station, it is critical to first characterize the wind regimes within and above a forest on the lee of a hill, and to determine the dynamical forces that control them. In this chapter, two series of numerical experiments and their results were discussed. The surface temperature during the day under a plant canopy is uniformly lower than over a bare soil surface. A dense canopy reduces the diurnal range of surface temperature because a large part of the incoming solar radiation is absorbed by plant surfaces near the treetop, reducing the amount of solar radiation reaching the ground and in turn decreasing heating to the neighboring air. A t night the outgoing long-wave radiation is also partly intercepted and absorbed by vegetation, which also radiates energy back to the surface. This  slows the radiative cooling of the soil. The combined effects of the two processes result in smaller diurnal range of temperatures on vegetative surfaces. The simulated wind flows at the West Coast Flux station show a distinct diurnal pattern above the canopy in fair weather. During the day the winds flow in the upslope direction, and at night the flows are mostly down-slope. But under the canopy, the statistics of simulated winds in the diurnal-cycle experiments indicate that down-slope flows occur 60% of time during 24 hours. The diurnal-cycle experiments captured the first-order effects of diurnal heating/cooling on the sloping terrain and its interaction with a heterogeneous rough surface. Above the canopy, the wind regime occurring at the West Coast Flux station is a combination of anabatic/katabatic flow plus the land/sea breezes. Inside the canopy the radiative effects of a tall dense canopy result in a relatively stable stratification near the surface during daytime. Meanwhile strong convective thermals in the upper canopy and vigorous turbulence mixing induced by the tree-drag were found to transport warm air from the canopy top. The combined effects of these factors result in intermittent down-slope flows even during the daytime, which mostly occur during the transition times and when strong convective thermals and turbulence dominate. The wind pattern near the sunset suggests a quick transition of wind direction and a reduction of wind speed. These phenomena result from a combination of the shadowing effects on the eastern slope (earlier sunset), and the west-slope dominates upslope winds extending over the ridge because of the stronger heating rates on the western slope. The rather rapid shift from up-slope to down-slope winds just after sunset may be an important factor in helping explain some of the CO2 flux anomalies just after sunset at the West Coast Flux tower. Compared to bare soil, the drainage flow for a dense canopy is much weaker and is separated into sub-canopy and above-canopy parts. The wind profiles indicate that the minimum wind speed occurs in the upper canopy, because most of the momentum is absorbed in the upper canopy where the leaves are concentrated. Farther down, as the leaf area density decreases to zero, tree drag gradually disappears, and the pressure-gradient term likely becomes more significant and dominates over the flux-divergence term. Consequently, the wind speed may increase with decreasing height above the ground, until surface friction reverses this trend close to ground level.  The experiments are useful for understanding the complex effects of tall trees on the airflow over a hill. The volume exchange of heat and momentum above and within a canopy located on a hill appear to be strongly affected by the local flows resulting from diurnal thermal forcing. The mechanisms causing the quick transition from anabatic to katabatic flow may help explain the anomalous flux measurements just after sunset.  Chapter 4 C 0 budget and advection effects 2  4.1 Introduction A better understanding of the global carbon balance, and particularly of the role of terrestrial ecosystems, presents an important challenge to the scientific community — the 'missing carbon sink' issue [Schimel, 1995; Ciais et al., 1995; Sarmiento and Grunber, 2002]. Given the importance of temperate forests in the global carbon balance [Reich and Bolstad, 2001], measuring and modeling the net-ecosystem exchange of CO2 (NEE) above these forests has become a major research activity. The most widely used method in observational studies of CO2 sources/sinks is the eddy covariance approach [Wofsy et al., 1993; Black et al., 1996, Saigusa et al., 1998; Holland et al., 1999; Carrara et al., 2003]. However, the eddy covariance technique has shortcomings under certain atmospheric conditions in which it does not represent the surface flux well. A t most sites, the sum of latent and sensible heat fluxes using the E C method is about 80-90% of the available energy flux even during the daytime [Aubinet et al. 2000]. This suggests that daytime CO2 fluxes may be underestimated. A t night, when wind speeds are often low, measured fluxes are less than respiration rates that are expected to occur. Even when the change in CO2 storage in the column of air beneath the E C sensors is taken into account, the flux is less than expected [Lavigne et al., 1997; Mahrt, 1998; Morgenstern et al., 2004]. This simply suggests that the respired CO2 must be transported away from its source by some other mechanisms that are undetected by the current method (storage + eddy flux) based on the assumptions of horizontal homogeneity and flat terrain. Many natural ecosystems are on complex terrain, and they usually form a non-uniform and patchy landscape pattern with discontinuities and border effects. The real conditions at most flux sites do not fully meet the requirements of E C approach. A few field experiments suggest that vertical advection could be significant over sloping terrain. For example, Eugster and Siegrist (2000) employed a tethered balloon to measure the vertical profiles of CO2 concentration and found that during daytime when the atmosphere was well mixed, the near-surface vertical gradient of CO2 was typically on the order of 2 to 5 ppm m" . A t dusk, 1  the vertical gradient of CO2 concentration decreased to around -20 to -25 ppm m" , and at 1  night it was about -10 to -20 ppm m" . So the vertical gradient of CO2 concentration always exists. If a non-zero mean vertical velocity exists, then vertical mean advection of CO2 should occur. Lee [1998] thought that a zero mean vertical velocity may be a good assumption very close to the ground, but it was questionable over tall vegetation. A t a height of 40 m above the ground (a typical height for observations) the mean vertical velocity may be significant. He suggested that neglecting vertical advection was the major error in the current N E E estimates, so he proposed a vertical advection correction to calculate N E E based on the scalar conservation budget in a one-dimensional framework. At the West Coast Flux station, the results in Chapter 3 suggest that three mechanisms could cause non-zero mean vertical velocity under fair weather. O f these, drainage flow on sloping terrain is the most important one for long-term flux observations, because the West Coast Flux station is located on a 5-10 degree slope. At this location, thermally driven mesoscale circulations such as anabatic/katabatic flows occur frequently, as shown in the observations and our simulations in Chapter 3. During fair weather daytime conditions, the convective turnover at a horizontal scale of 2 ~ 4 km could also contribute to the non-zero vertical advection [Garratt, 1992]. In such situations, vertical velocities as large as 10 cm/s are possible at the height of the tower (42 m above the ground) within the descending air column, according to our simulations and the study by Moeng [1984]. Synoptic-scale subsidence could be another mechanism causing non-zero mean vertical velocity under fair weather. For a heterogeneous surface, mean horizontal advection could also be significant. To date, only a few field experiments have addressed horizontal transport of CO2 with limited data, and no published information is available on numerically modeling advection effects on flux measurements and N E E estimates. Sun et al. [1998] indirectly inferred the significance of horizontal advection by showing increased CO2 concentrations over a lake at night and in the early morning due to nocturnal advection of respired CO2 from the surrounding forests by the land-breeze circulation. The horizontal transport was expected because of the large horizontal gradient of CO2 concentration for spatially heterogeneous vegetation — from forest to lake. In a study site located in a heavily forested zone with a very tall tower located in a grassy clearing of -180 m radius, Y i et al. [2000] found  heterogeneous vegetation leading to advection accounting for 27% and 5% of diurnally integrated N E E values between levels of 30 and 122 m, and 122 and 396 m, respectively. Based on tethered-balloon measurements from two nights, Eugster and Siegrist [2000] attributed the difference between the fluxes estimated by the boundary-layer budget method and the surface eddy fluxes to be horizontal advection. They found the advection was 12.5 times that of the local vertical eddy flux. Staebler and Fitzjarrald [2004] found that including horizontal transport didn't fully account for the observed difference in N E E between flux deficit and non-deficit nights, but decreased the difference. For the West Coast Flux station, Georgia Strait is only 9 km away, providing a welldefined surface heterogeneity and preferred CO2 venting location in the nocturnal A B L . Under clear nocturnal skies, radiative cooling of the land surface and nearly timeindependent temperature of the water surface leads to differential surface cooling and generation of the land breeze, as shown in Chapter 3. When the synoptic wind is weak, the land breeze and drainage flow could effectively redistribute CO2 locally from the land to the water surface. If CO2 is transported away horizontally or vertically and is undetected by the E C technique, then the carbon budget could be significantly affected. Furthermore, as we mentioned before, most other E C sites have also reported this problem, i.e., reduced turbulent mixing at night and underestimation of N E E . Therefore, the real condition at most E C sites are not as ideal as required by the E C approach, and the complex topography at the tower site raises uncertainties on E C measurements and N E E estimates. In this chapter, we focus on simulating the diurnal response of CO2 field and the advection effects on E C measurements, and analyzing the contribution of each term in the CO2 budget and the N E E estimate. The incorporation of the CO2 budget equation into WFIS with a treatment of sources/sinks is introduced in Section 4.2, followed by a description of two numerical experiments, results and discussion in Section 4.3. Using the outputs of all meteorological variables and the mixing ratio of CO2 from the numerical experiments, the contributions of all terms in the CO2 budget are quantified and advection effects are assessed in Section 4.4. Conclusions are given in Section 4.5.  4.2 Treatment of C 0 budget in the model 2  The conservation equation of C 0 provides the foundation for calculating CO2 fluxes: 2  (4-1) where c represents the mixing ratio (kg/kg) of carbon dioxide, p is the air density (kg/m ), a  u, v and w are the wind velocity components along x, y and z directions, respectively. The over-bar in (4-1) denotes the grid (volume) average in the model, and is different from the time average used in the eddy-covariance measurements. The primes in (4-1) denote the sub-grid scale terms. If we only consider two dimensions and neglect the molecular diffusion term (D ) at the grid scale, (4-1) can be simplified to: (4-2)  •) = S (y,z) c  (VI) Equation (4-2) states that the local time rate of change of CO2 in a control volume (I) is balanced by the mean horizontal advection (II), the horizontal turbulent flux divergence (IV), the mean vertical advection (III), the vertical turbulent flux divergence (V) and the biological source/sink (VI). The turbulent flux terms (IV) and (V) in a numerical model essentially are treated as a sub-grid scale term and are parameterized with different closure schemes. To incorporate the CO2 budget (4-2) into the model, the most challenging problem is how to parameterize the biological source/sink (VI) at the various vertical levels in the model. The basic source and sink processes of CO2 occurring in an ecosystem are shown in Fig. 4.1. The major sink process of CO2 in a forest is the canopy gross photosynthesis (P), which consumes CO2 from the atmosphere and results in carbon storage in leaves, stems and roots. The source processes of CO2 include the autotrophic respiration (R ) by leaves, a  boles and roots, and the heterotrophic respiration (R ) from the litter-fall and the soil, which s  return CO2 to the atmosphere. The sum of the sources/sinks of CO2 from an ecosystem to the atmosphere (NEE) is given as: NEE = R + R a  s  (4-3)  -P  It is generally accepted that the soil CO2 efflux (R ) often accounts for more than 70% s  of forest ecosystem respiration [Janssens et al., 2001; Law et a l , 1999]. Considering the uncertainties of the autotrophic respiration (R ) from the boles, leaves and roots, we assume a  all the ecosystem respiration occurs right above the ground and ignore the autotrophic respiration (R ) by leaves and boles, so that a  NEE = R -P  (4-4)  S  Canopy Carbon Balance Net Ecosystem Carbon Exchange  U ^ till)  Gross Photosynthesis  ;mti  Dark and ^fc? J Photo Respiration • m Xp^^CO^Storage  I  lobule respiration  root respiration  microbial respiration  li'tlei  -  ™M,iiraHi.n  Figure 4.1: Schematic diagram of carbon sources/sinks http://namre.Berkeley.Edu/biometlab].  in  an  ecosystem  [from  It has been found that the CO2 flux from the forest floor is highly correlated to the soil temperature [Law et al., 1999; Janssens et al., 2001; Drewitt et al., 2002; Jassal et al., 2004]. To parameterize ecosystem respiration, the most common choice is to use a functional relation between soil temperature and soil C O 2 efflux, which is often an exponential function. In the model, the empirical formula obtained from the soil chamber measurements at the West Coast Flux station is used to quantify ecosystem respiration R as s  Drewitt [2002] and Morgenstern [2004]: R = R Q^ "' T  s  where  i? f is the respiration rate at a reference re  (4-5)  )la  re/  temperature T f. In the model, a - 10 °C, K  rf re  = 10 °C and R { =5.0 ^mol m"V as suggested by Morgenstern [2004]. Q\o is the factor by Te  which R increases for a 10 °C increase in soil temperature [Lloyd and Taylor, 1994]. The s  value of Q\o is in the range of 2~3 during the warm season (in July) as reported by Drewitt et al. [2002], we use Q  l0  = 2.2. The soil temperature at 5 cm depth is obtained using a 4-  layer soil heat transfer model in the WFIS. The respiration calculated by (4-5) has flux units: u.molm'V . 1  The most widely used relationship for describing the light response of photosynthesis (P, fimol m ' V ) at the canopy scale is the Michaelis-Menten function [e.g., Chen et al, 1  1999; Morgenstern et al, 2004]: p = "AyooQ where a is the quantum yield  (4-6)  /jmo\ CO (a value of 0.1 was used in the model), A fjmol photons 2  represents the canopy photosynthetic capacity [fimol m " V ] , a value of A = 1  MM  V  1  R  22 |amol m"  was used, Q is the down-welling photosynthetic active radiation (PAR) [ ^ ° ^ ^  o X o n s  ^  ms which is the visible part of the short wave (SW) radiation (0.38 - 0.71 u.m). It is directly related to the incoming SW at the top of the canopy, and can be roughly calculated as two times SW. The budget of incoming solar radiation and net radiation at the canopy top for a four-day simulation NOF2 experiment described in Chapter 3 are shown in Fig. 4.2.  -—  5  10  20  30  AO  SO GO T i m e (Hours)  solar  7Q  SO  90  10O  Figure 4.2: Radiation budget at the canopy top for a 4-day simulation NOF2.  Since (4-6) gives the whole ecosystem photosynthesis, we have to parameterize it into multi-layers in the model using leaf area index similar to the C A N V E G model [Baldocchi, 1998]. The parameterization scheme for photosynthesis in each vertical layer down the whole ecosystem photosynthesis of (4-6) as  is to scale  = -/>{-!- 1 a(z)dz} * S&L  P(z ) k  LAI o  (4-7)  C{z _{) k  where P (umol m'V) is calculated by (4-6), a(z) (m /m ) represents the leaf area density. 2  3  LAI is the leaf area index, C(zk) (kg/kg) is the CO2 mixing ratio at height z (at grid level k), and C(zk-i) is the CO2 mixing ratio at the lower grid level (k-\). The reason for multiplying the ratio (C(zk)/C(zk-i)) is to avoid a negative mixing ratio of CO2 being generated in the model. Photosynthesis is an atmospheric sink for CO2, and i f we simply split (4-6) into layers and apply it, we might end up with negative CO2 mixing ratios. Substituting (4-6) into (4-7) and (4-5) into (4-4), we can get the net ecosystem source flux of CO2 (NEE). The source/sink in (4-2) is the gradient of N E E : s  c  ( y , ^ ^  (4-8)  az  where N E E has units: |iunol m ' V , and the gradient of N E E has units of umol m ' V . Thus 1  1  the units on both sides of the CO2 budget equation (4-2) match.  4.3 C 0 budget experiments and results 2  Two experiments were performed to investigate the diurnal response of CO2 fields to localscale atmospheric motions and to test the code of CO2 budget in the model. Table 4.1 outlines the characteristics of the CO2 budget experiments. A l l symbols have the same meaning as described in Chapter 3. These experiments adopt the same domain configuration, topography, and canopy specification as those in the diurnal cycle experiments (Fig. 3.1). The dynamics and physics in the model were all active with the CO2 budget active. The initial CO2 field was set to 370 ppm (0.561 g C C V k g air). We used periodic boundary condition for CO2 field. The simulation duration for the two experiments was 30 hours starting at 0500 LST. Zero background wind and clear sky were assumed.  Table 4.1: Characteristics of C 0 budget experiments. Ay and Az represent the grid size of horizontal and vertical, respectively. At denotes the time step. Slashes separate values for the two nested domains. Exp. Az \ Ay(m) H i l l height Duration At (s) C0 (m) (m) budget (hours) (tree-cover) CO2F0 150/75 12/6 0 1300 1.0/0.5 On 30 2  m n  C02F1  150/75  12/6  2  3.0/1.5  1.0  1300  On  30  4.3.1 Simulated soil temperatures Time series are shown in Fig. 4.3 for the simulated soil temperatures at several depths for the dense canopy case (C02F1, r\ = 1) and bare soil case (CO2F0, r\ = 0). The evolution of soil temperature is for grid point 392 (located at y = 58.8 km over the land) in the outer domain, which has the same elevation (320 m) as the West Coast Flux station. For bare ground (r) = 0), the range of soil temperature is much larger than that for denser canopy. A s we learnt in Chapter 3, a dense canopy reduces the diurnal range of soil temperature because a large part of the incoming solar radiation is absorbed. Hence, during the day the soil temperature under a plant canopy is uniformly lower than that over a bare surface, and is warmer at night. Due to the rapid radiative cooling at night, the soil temperatures in the upper layers are cooler than those of the deep soil layers for the bare-soil case. Under a dense canopy, the outgoing long-wave radiation is partly intercepted and absorbed by vegetation, which also radiates energy back to the surface, this combination of effects slows the radiative cooling of the soil. The soil temperature in the upper layer is a little warmer than that of the deep soil at night. The soil temperature is almost constant with time at a depth of 59 cm. 296  4.3.2 Simulated CO2 flux by respiration and photosynthesis Substituting the simulated soil temperature at the 5 cm depth into the empirical formula (45), the simulated CO2 flux by soil respiration for experiments C 0 2 F 1 and CO2F0 are shown in Fig. 4.4. Higher soil temperatures result in stronger soil CO2 efflux. The simulated CO2 fluxes by canopy photosynthesis, using empirical formula (4-6) and parameterization scheme (4-7), are shown in Fig. 4.5. The difference is caused by multiplying the ratio (C(z )/C(z .,)) in (4-7). k  k  20 T i m e (hrs)  25  Figure 4.4: Simulated soil respiration for experiment C02F1 (upper panel) and CO2F0 (lower panel).  simulated using parameterization scheme in the mode! - — calculated by empirical formula  .S3  •°  E  I .1 o = o  kti 20  T i m e (hrs)  Figure 4.5: Simulated photosynthetic C 0 flux using empirical formula (4-6) (dash line) and parameterized scheme (4-7) (solid line). 2  4.3.3 Simulated CO2 profiles in the A B L Six simulated vertical profiles of CO2 concentration at grid point 392 (y=58.8 km on the land) of the outer domain show a diurnal sequence of the CO2 field (Fig. 4.6). Grid location 392 in the outer domain corresponds to the location of grid point 280 in the inner domain. During daytime there is a well-mixed profile from the top of the A B L down to the canopy top due to strong convection and turbulent mixing (1100, 1500 and 1900 L S T in Fig. 4.6). A minimum CO2 concentration over the canopy top and a steady decrease of CO2 concentrations in the lowest 200 m of the A B L from 0700 to 1900 L S T is a result of ongoing photosynthetic uptake of CO2 by the canopy. The minimum CO2 concentration during the day is found at about the 400 m height just before sunset (1900 LST) when the energy balance at the canopy top begins to change from positive to negative; i.e., the net long-wave radiation loss exceeds the short-wave radiation gain at the canopy top. Around sunset, CO2 concentrations near the surface start to build up in a very thin layer, resulting in a very strong negative vertical gradient of CO2 near the surface. As time progresses the CC^-rich layer at the surface deepens and the CO2 concentration reaches its maximum value, around 650 ppm near the surface at midnight (2300 LST), the depth of the C O v r i c h layer reaches to roughly 400 m above the ground at 0300 L S T . From midnight to sunrise, there is no increase in CO2 concentration at the ground, and it even decreases from 650 ppm to 620 ppm. Based on observation data, Yang et al. [1999] also reported the phenomena of decreasing tendency of CO2 concentration after midnight. If the CO2 originates from only local soil respiration, a steady increase in CO2 concentration should be expected. Then where is the respired CO2 going? W e ' l l discuss this later. At dawn (0700 LST), a distinct change in the vertical distribution of CO2 concentration reveals the height of the stable nocturnal boundary layer (« 150 m above ground), which is now affected by the exchange processes that are switched on when the solar radiation quickly changes the surface energy balance. The profile at 0700 L S T shows two peaks: one minor peak at a height of 150 m indicates warmer C02-rich air is being advected over the cold layer near the surface by the onset of upslope flow. Another CO2 peak near the surface is decreasing from 600 ppm to 450 ppm, which is a combined effect of the onset of canopy photosynthesis and the beginning of atmospheric mixing after sunrise. After a transition  period, strong mixing eliminates the steep gradient near the surface, and the CO2 profile becomes vertically uniform (1100 and 1500 LST).  1000  1100  1000  800 600  f  o  400  o 200  1000  800  800  600  600  400  400  200  200  0 300 400 500 600 700  0 300 400 500 600 700  1000  1500  1000  2300  1900  0  = • © - — 0 300 400 500 600 700  0300  800 , o  _ 800 E  600  •£ 600  o  400 ® 400  () 200  200 0300 400 500 600 700  0 300 400 500 600 700  C0  2  300 400 500 600 700  concentration (ppm)  Figure 4.6: Simulated vertical profiles of C 0 concentration at 6 times (LST) during a day at grid location 392 (at 58.8 km on the land). 2  To address the decreasing tendency of CO2 concentration on the land after midnight (Fig. 4.6) and assess where the respired CO2 goes, we chose another grid point 450 at location 67.5 km of the domain to see the time evolution of CO2 vertical profiles (Fig. 4.7). This grid point is located over the water of Georgia Strait. During daytime, the CO2 field is well mixed (1100 LST). The CO2 concentrations are nearly constant down to the water surface. However, in the upper boundary layer (8001000 m), a weak decrease of CO2 occurs because the land breeze (the return flow of the  sea breeze) in the upper A B L is advecting the COvpoor air from Vancouver Island to Georgia Strait (1500 LST). No minimum CO2 concentration can be seen near the water surface because there is no canopy consumption. As time progresses,  the C O 2  concentrations also start to accumulate near the water surface after sunset, and the accumulation process lasts throughout the night. Until 0300 LST, CO2 concentration over the water surface has a similar magnitude as that on the land at this time. Since no other mechanisms locally emit CO2 over Georgia Strait in the model, it can only be transported from the neighboring land surface.  1000 800  1000  1000  1  1100  i  |  800  800  •3i  <=- 600 -4—'  si  a> CD 400 I 200  1500  ir 'r  t  1900  t I  600  600  •-  I  400  i t i  1t  400  •S•-  i  200  i i  i.  200  . . ,^  t e •  t. . . 300 400 500 600 700  i. . . 300 400 500 600 700  0 •n— 300 400 500 600 700  300 400 500 600 700  300 400 500 600 700  300 400 500 600 700  C 0 concentration (ppm) 2  Figure 4.7: Similar to Fig. 4.6 except for grid point 450 (at location 67.5 km over the water surface).  4.3.4 Diurnal variations of CO2 concentration and gradients Figure 4.8 shows the diurnal variation of CO2 concentration and its vertical gradients for the lowest three layers (6 m, 18 m and 33 m) for experiments C 0 2 F 1 . Figure 4.9 is similar to Fig. 4.8 except showing CO2 concentration and its gradient at 3 m, 9 m and 15 m above ground for the inner domain. During daytime when the atmosphere is convectively unstable, the CO2 concentration remains nearly steady with time due to strong eddy mixing, while the ongoing photosynthetic uptake of CO2 by the dense canopy results in a very slight but steady decrease (in upper panels of Fig. 4.8 and 4.9). At night CO2 accumulates near the surface because the strong stability above the dense canopy suppresses upward turbulent transport of CO2, and traps under the canopy the respired CO2 by soil. During daytime, the vertical CO2 gradient near the surface is typically on the order of -2 ppm m" between 6 m and 18 m above the ground (lower panel of Fig. 4.8), and -5 ppm m" 1  1  between 3 m and 9 m (lower panel of Fig. 4.9). The gradient across the 18-33 m layer is much smaller and close to zero or even positive due to the combination of photosynthesis and strong mixing induced by the dense canopy. Around sunset, when the atmospheric stability rapidly increases, the vertical CO2 gradient increases quickly for a one-hour period (we call it a transition time). The average value for the period is around -10 —15 ppm m"  1  between 6 m and 18 m (lower panel of Fig. 4.8); and around -20 to -25 ppm m" between 3 1  m and 9 m (lower panel of Fig. 4.9) and is in good agreement with the observations [Eugster, 2000]. The rapid increase in CO2 concentration ceases at about 2000 L S T but continues to slowly increase to midnight. After that, the high CO2 concentration remains relatively constant for the rest of the night, indicating a quasi steady state with the rate of losing CO2 (through advection processes) being approximately equal to the efflux of respiratory CO2 from the soil. After the transition period, the vertical CO2 gradient quickly drops back to values around -5 ppm m" in the 6-18 m layer and around -10 ppm m" in the 3-9 m layer. 1  1  Another peak of CO2 gradient shows up around sunrise corresponding to the large drop in near-surface CO2 concentration in the early morning. This phenomenon of CO2 rapidly dropping in the early morning is also commonly observed in dense forests [Grace et al., 1995; Yang et al., 1999].  The near-surface vertical CO2 gradient slowly but steadily decreases throughout the whole night (arrows in lower panel in Fig. 4.8 and 4.9), although the high C 0  2  concentrations persist during the period. Both a steady (or slightly increasing) C O 2 concentration and gradient should be expected i f the C O 2 originates only from local respiratory processes, suggesting there must be other processes  such as advection  accounting for the decreased CO2 gradient. Further, eddy covariance flux measurements (Fig. 1.3) at the West Coast Flux station show a significant peak of upward eddy flux of CO2 around sunset, which coincides with the large C 0 gradient occurring at this time in the simulations (lower panel of Fig. 4.8 and 2  4.9). The simulated C 0 concentration profiles (1900 L S T in Fig. 4.6) capture the correct 2  time of change from the daytime with C 0 uptake by the forest to the nocturnal regime with 2  CO2 release from the soil. At this time, a steep gradient of CO2 exists only in the near surface (1900 L S T in Fig. 4.6), while above the canopy the CO2 gradient is close to zero. For the bare-ground case (Fig. 4.10), the C O 2 concentration remains  constant  throughout the day due to the stronger thermal convection and turbulence mixing under fair-weather conditions. At night the C 0  2  concentration accumulates near the surface  similar to the dense-canopy case, but has a smaller magnitude — less than 500 ppm. Two processes could account for the difference. One is the decreased soil respiration. Since the soil respiration is mainly a function of soil temperature, radiation cooling is much stronger over bare ground and may result in a lower soil temperature and hence a lower soil CO2 efflux. Another is much stronger advection by cold drainage flow over bare ground (Fig. 3.7) that transports the respired CO2 away from where it originated. The CO2 gradient for the bare ground case also shows peaks around sunrise and sunset and is similar to the dense-canopy case. Other than at transition time, the vertical C 0  2  gradient during day is similar to that at night for the bare-ground case, with values around 4 ppm m" in the 6-18 m layer and -0.5 ppm m" in the 18-33 m layer. These values are 1  1  much smaller than those under the dense-canopy because of the stronger daytime thermal convection and nighttime drainage flow over bare ground (Fig. 3.7). The time series of horizontal C 0 gradients show intermittent CO2 variation from -10 to 2  10 ppm per 150 m during daytime, and steady positive values of 2 ppm per 150 m at night (Fig. 4.11). This indicates that although we assume a horizontally homogenous canopy  covering the land in the model, the sloping topography (varying elevation) results in varying soil temperatures that lead to varying CO2 respiration (a value of 2 ppm per 150 m). During the day the strong convective thermals increase the variation.  Local time (hrs) Figure 4.8: The time series of C 0 concentration (upper panel) at the height of 6 m, 18 m and 33 m, and C 0 concentration gradient between the two layers (lower panel) for the outer domain of the dense-canopy experiment C02F1. Nighttime is during about 19 to 29 hrs local time. The arrow indicates the decreasing tendency of C 0 gradient between 6 m and 18 m at night. 2  2  2  p 800  o o 200 • CM  0 0  o  1  5  >-  1 15  10  1 20  1 25  1  1 30  35  Local time (hrs) Figure 4.9: Diumal courses of C 0 concentration (upper panel) at the height of 3 m, 9 m and 15 m, and CO2 gradient between these two vertical layers (lower panel) for the nested domain of the dense-canopy experiment C02F1. The arrow indicates the decreasing tendency of C 0 gradient between 3 m and 9 m at night. 2  2  800 c o  700 -  bare ground  ra  600 c "E 500 a> o a , c o — 400 o CN 300 O o 200 i—  — 6m — 18m - 33m  100 0  10  15  20  25  Time (hrs)  30  35  30  35  -15  +5  — 10  gradient between 6m and 18m gradient between 18m and 33m 15  20  25  Time (hrs) Figure 4.10: Same as Fig.4.8 except for the bare-ground case (CO2F0).  15  T i m e (hrs)  20  25  29  Figure 4.11: Time series of horizontal C 0 gradient at 3 m high above ground at grid point 280 for the dense-canopy case. 2  4.3.5 CCh spatial distribution The simulated CO2 spatial fields are shown in Fig.4.12 and display advective processes superimposed with the photosynthetic and respiratory processes of the ecosystem. We start at midday (1300 LST) with a well-mixed atmospheric condition. The C 0 concentration 2  near the surface does not differ much from that aloft due to the strong turbulent mixing in the convectively unstable A B L . The early sunset on the east side of the mountain at 1900 L S T turns the surface energy balance negative earlier than on the western slope. The thermal convection ceases and down-slope flow of C02-rich air first begins to form on this side (Fig. 3.5). Meanwhile the west side of the mountain remains under sunlight, the upslope flow dominates and the CO2 field is still well mixed. Thermal convection ceases on both sides of the mountain at 2100 L S T , and the CO2 concentrations build locally on the land surface. Around midnight (2300 and 0100 L S T ) a near-surface accumulation of CO2 concentration has formed both on the land and over the neighboring water surface. This is an indication that C 0 - r i c h air is being transported down 2  to the neighboring water surface where CO2 concentration is initially much lower. Under clear skies with weak synoptic flow, the results in Chapter 3 (Figs. 3.5, 3.6, and 3.7) showed that the dominant flow at night is drainage flow. The flow is further amplified by the nocturnal land breeze at the West Coast Flux station (Fig 3.8), because of the heterogeneous surface between Vancouver Island and Georgia Strait. It is the drainage flow that is advecting the CCVrich air from the land to the water surface and causing the CO2 to accumulate over the water. The accumulation of CO2 over the water surface reaches its maximum depth before dawn (0300 and 0500 LST). The early sunrise on the east side of the mountain first breaks up the stable boundary layer, and the accumulated C 0 concentration flushes out (0700 2  LST), leading to a rapid fall of CO2 concentrations near the surface. The time evolution of CO2 spatial distribution (Fig. 4.12) shows the capability of nocturnal drainage flow and weak land breezes to advect the locally respired CO2 on a regional scale. The underestimation of nocturnal respiration based on eddy covariance measurements is likely associated with the advection by drainage flow. After sunset, the thermal convection ceases resulting in a calm evening with weak turbulence. Advected C 0 by the drainage flow becomes a dominant portion of the CO2 budget.  2  17.3  34.6 Y  51.8  69.1  86.4  I KM)  Figure 4.12: Simulated C 0 spatial distributions at different times in a day for dense-canopy experiment C02F1. 2  In summary, the simulated vertical profiles and time series of CO2 concentration show a well-mixed CO2 field with a minimum CO2 concentration in the upper canopy during the day. A rapid build up occurs in the trunk space around sunset. The C O v r i c h layer deepens and reaches its maximum depth near midnight. For the rest of the night the C O 2 concentration slowly decrease. The near-surface CO2 accumulation over the neighboring water results from CO2 advection from the land at night. The vertical CO2 gradients above 9 m are small. Below the 9-m height there are significant vertical gradients and diurnal variations, particularly around sunset and sunrise. The simulated horizontal CO2 gradients are significantly smaller than the vertical gradients. The simulated spatial distributions of CO2 field at various times (Fig. 4.12) are used to interpret the various processes that are active. At night the drainage flow and weak land breeze advect the locally respired CO2 from the land to the water surface. This process explains the slowly decreasing C O 2 concentration and its gradients near the surface after midnight.  4.4 C 0 budget analysis 2  In this section we utilize the meteorological field variables (v, w, K ) and CO2 (c) from m  experiment C 0 2 F 1 to quantify the contributions of each term in the CO2 budget. The N E E is calculated from the integration of each term in the CO2 budget from the ground to the tower height. The variables such as the wind velocities (v and w), CO2 mixing ratio (c) and eddy coefficient (K ) are output every 15 seconds from the numerical simulations. m  The WFIS model treats the sub-grid-scale (SGS) flux term in (4-2) as p W?=-p a  a  K dc Idz m  (4-9)  which is a typical first-order closure. A similar approximation is also used in the y direction. K is eddy-mixing coefficient and calculated in the model as Smagorinsky [1963] and L i l l y m  [1962]. The over-bar represents a grid-volume average in the model, and is different than the time-average used in the E C technique. These differences need to be kept in mind when comparing the resolved-scale vertical fluxes from the model with the eddy-flux measurements at the tower site. How one computes the mean and fluctuations from the mean is an issue that affects calculations using equation (4-2). The assessment of the flux covariance requires we sample all significant scales of motions contributing to the turbulent transfer of CO2 in the  atmosphere. The current sampling rates used in the E C technique (10 to 20 Hz) might enable one to capture the high-frequency portion of the flux cospectrum, but to capture lowfrequency contribution associated with the daytime convective boundary layer and nocturnal stable boundary layer, the practical time average of turbulence over a 30 minutes period for day and night is questionable [Lee, 1996; Massman and Lee, 2002]. It is still an unresolved issue that how long a time average should be used for the eddy-flux measurements in different atmospheric conditions [Baldocchi, 2003]. Thus, numerical studies could be used to assess optimal averaging techniques to further improve the E C method The sum of the mean and turbulent fluctuation terms in the model constitutes the estimates of both the horizontal and vertical advection. The vertical CO2 flux is composed of the two terms: QQ-  F = p w"c+~pw'c' = p Wc - p K — a  (4-10)  m  dz  which denote the resolved vertical fluxes and sub-grid-scale (SGS) contribution. Using horizontal and vertical integration, one can formulate various volume tendencies for the atmospheric storage of CO2 along with the contributions from horizontal and vertical advection. Substituting (4-9) into (4-2), we get: d(p c)  d(p  —^-+[—^ n  vc)  dt  dy  (I)  (II)  d(p  -+_^  Wc) dz  d(-p  -] + [  (III)  P  K  a  m  dcldy)  dy  d(-p  —+  a  K m dc I dz)  e  dz  (IV)  _  -] = S (y,z)  (V)  (4-11)  (VI)  Thus, N E E is essentially an integration of each term on the left side of (4-11) from the ground to the height of the tower instruments (42 m): hd{p-c)  NEE = [ f o  t>d(pjc)  dz + f  a  a  [1]  Q  d  hd(pWc)  y  [2]  hd(-pK dc/dy) m  dz +  a  0  d  +I z  [3]  Q  d  a  m  h d(-p  -dz] + I 0  y  [4]  K dc I dz) a  m  -dz  ,  A  (4-12)  &  [5]  4.4.1 Contributions of each term in the CO2 budget Figure 4.13 shows the contribution of each term on the left side of CO2 budget (4-11) at 1300, 1900, and 0100 L S T (half-hour average), respectively. At 1300 L S T the storage term (I) is very close to zero due to the strong eddy mixing. The positive values of the resolved  horizontal advection (term II) indicate CO2 divergence causing CO2 loss in a volume. The resolved vertical advection (term III) has negative values above 6 m denoting CO2 is mixing into the volume from upper canopy to compensate for the C O 2 consumed by canopy photosynthesis. The contributions from the horizontal SGS flux gradient of CO2 (IV) are very small compared to other terms. Thus, it is safe to neglect this term. The vertical SGS flux term is fluctuated at different height. The resolved vertical and horizontal advections have opposite signs and are the biggest terms at this time. We can't ignore any of them to correctly estimate the source. At 1900 LST, the contribution of each term is considerably smaller than that at 1300 L S T . The storage term shows positive values in the trunk space and negative values in the upper canopy, indicating that CO2 starts to build up in the trunk space while canopy photosynthesis is still actively consuming CO2 in the upper canopy. Negative horizontal advection (term II) is bringing CO2 into the volume while positive vertical advection (term III) is transporting the CO2 away. The horizontal SGS flux (term IV) is still very small. The vertical SGS flux (term V ) is showing CO2 loss in the upper canopy and CO2 gains from the bottom of the volume. The different sign of the storage term in the upper canopy to that in the trunk space indicate that this is a transition time. A l l terms are one magnitude smaller near midnight (0100 L S T ) compared to daytime, due to a calm and stable nighttime A B L . The storage term has small negative values at all vertical levels indicating a slight loss of CO2 from the volume at this time. The resolved horizontal advective fluxes (term II) and the vertical SGS fluxes (term V ) are responsible for transporting CO2 away from the volume, while resolved vertical fluxes (term III) is bringing CO2 in. To see the combined effects of the resolved fluxes and sub-grid-scale fluxes both in the horizontal and vertical directions, Fig. 4.14 shows the contributions of the storage (term I), horizontal advection (term II+IV) and vertical advection (III +V) at 1300, 1900 and 0100 L S T , respectively. A t 1300 L S T the horizontal advection is responsible for decreasing CO2 in the volume while vertical advection accounts for a CO2 gain; the net effect is near zero storage of CO2 due to strong convection and eddy mixing. A t 0100 L S T the horizontal advection is transporting CO2 away, and vertical advection below 6 m is also advecting  CO2 away; but above 6 m vertical advection accounts for the accumulation of CO2, the net effect is a slight loss of CO2 after midnight.  40 -e- storage I - B - horizontal advection II -0- vertical advection III horizontal SGS advection IV vertical SGS advection V  30  s  D)20 Z  10 0100 LST -0.2  0.2  0.4  0.8  0.6  1.2  C 0 flux gradient (micromole/m /s) 2  Figure 4.13: Simulated contribution of each term in the C 0 budget equation (4-11) at 1300, 1900 and 0100 LST. The symbols are circles (storage term I), squares (mean horizontal advection term II), diamonds (vertical advection term III), dot line (horizontal SGS flux term IV), and asterisks (vertical SGS flux term V). 2  Around sunset (1900 LST)  the situation is complex. In the upper canopy all terms are  negative indicating a C O 2 gain. Below 20 m, the horizontal and vertical fluxes have opposite signs and largely counteract each other. The CO2  storage term is negative at the  top of the volume and positive in the bottom, indicating that CO2 space while canopy photosynthesis is still actively consuming CO2  -0.2  -0.1 0  builds up in the trunk in the upper canopy.  0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8  CO flux gradient (micromole/m /s) 3  Figure 4.14: The contributions of horizontal advection (squares, term II+IV), vertical advection (asterisks, term III+V) and storage (circles, term I) in C02 budget (4-11) at 1300, 1900, and 0100 LST.  - 4 - 3 - 2 - 1  -0.4  -0.2  0  0  0.2  1  0.4  2  3  0.6  0.8  CO flux gradient (micromole/m /s) 3  Figure 4.15: Estimation of the net source/sink using the sum of different terms. The symbols are asterisks (source estimated by the sum of horizontal flux gradient and storage), diamonds (the sum of vertical flux gradient and storage), circles (the sum of five terms), dot line represents the real source we incorporated in the model. Figure 4.15 shows the estimations of the net source/sink on the right side of the  CO2  budget (term VI) using the sum of different terms on the left side of (4-11). The real source term (dot line) is computed using empirical formulas (4-5) and (4-7) [Morgenstern, 2004]. It is evident that the sum of five terms is the best representative of the real source as the horizontal and vertical advections generally have opposite signs and considerably counteract each other. The sum of three terms is the worst because the mean horizontal  advection is needed to balance the mean vertical advection. At daytime (1300 L S T ) and nighttime (0100 LST), the sums of three terms all underestimate the soil respiration (source) near the surface, around sunset (1900 LST), the sum of storage and vertical flux underestimates the sink and overestimates the source (soil respiration).  4.4.2 Contributions of each term in N E E estimation One obtains the contribution of each term in the N E E budget (4-12) by integrating each term in CO2 budget (4-11) from the ground to the height of the flux instruments (42 m). A half-hour time filter (120-point running average) is performed for the model data to filter the noise. The average diurnal behavior of these terms and their variability are illustrated in Fig. 4.16. The storage term [1] shows a variation between -1 to +1 (j.mol m'V during the day 1  and a build-up of CO2 in the first half of evening (positive storage). A particularly rapid accumulation occurs around sunset when strong thermal convection ceases. Near zero or slight negative value of CO2 storage after midnight indicate a relatively steady state with the rate of CO2 loss approximately equal to the respired CO2 from the soil. A rapid exhalation happens around sunrise when the stable boundary layer breaks up. The resolved horizontal (term [2]) and vertical advection (term [3]) are negatively correlated. The advective fluxes are an order of magnitude larger than storage term during the day and display their largest variation during active convection. At night when thermal convection ceases, the resolved vertical advection decreases to near -1 (amol m'V . The 1  resolved horizontal advection has a positive magnitude around 2 - 4 (imol m"V, which is responsible for transporting locally respired CO2 away. The contribution of horizontal SGS flux (term [4]) is much smaller than other terms and plays a minor role in the N E E budget. The vertical SGS flux is negative during daytime and moves CO2 into the volume to compensate the CO2 loss by the canopy photosynthesis. Around sunset when the vertical gradient of CO2 is the largest, the vertical SGS flux grows rapidly to a magnitude around 5 ~ 6 umol m"V. For the rest of the night, it is the dominant term in the vertical direction. The positive values indicate that CO2 is vented away.  In the model, the sum of the resolved and SGS fluxes constitute the estimates of the horizontal and vertical advection (4-10), which are shown in Fig. 4.17. It is evident that during the day they have large amplitudes and opposite signs, and mostly counteract each other. However, it is interesting to see they have the same sign at night, i.e. positive values at night (after 2100 LST), which indicates that CO2 is advected away both by horizontal advection and vertical advection at night. The magnitude of advected CO2 flux is almost equal to the locally respired CO2 (Fig. 4.2) by the soil after midnight. That explains why the CO2 concentration doesn't increase during the second half of the night. The results of using the sum of all five terms and the sum of three terms (storage + vertical advection) on the right side of (4-12) to estimate N E E are shown in Fig. 4.18. During the day the ecosystem is clearly a sink, and at night it is a source. The time integrations would be a CO2 sink because the negative area clearly exceeds the positive area in Fig. 4.18. We compared the incorporated N E E using the empirical formula [NEE = respiration (^-photosynthesis (P), which is the source we incorporated into the model] to the estimated N E E with the sum of five terms (4-12). They are similar (Fig. 4.18a), indicating that the sum of all five terms is the best representation of the source. If we only use the sum of storage, resolved vertical advection and vertical SGS fluxes [terms [1] + [3] + [5] in (4-12)] to estimate N E E (Fig. 4.18b), it is evident that it underestimates by about 40% of the real source at night, and has large fluctuation during the day due to the strong convection and the lack of horizontal advection to offset the vertical advection. The E C technique samples turbulent motions using fast-response instruments (20 HZ). The vertical eddy fluxes are estimated by statistical analysis of these instantaneous vertical flux measurements using half-hour time averaging [Baldocchi, 2003]. To mimic this procedure currently employed by Fluxnet researchers to calculate vertical eddy flux, we perform a similar analysis procedure on the model output data (vertical velocity w and CO2 mixing ratio c every 15 seconds). We first average field variables (w and c) every 30 minutes to obtain w and c , and then the fluctuations (w' and c') are obtained individually (as is done in the E C technique) by subtracting the half-hour means from the data; i.e., W = w- w, c' = c-c  . The resolved vertical eddy fluxes are calculated by p w'c'. Then the a  N E E is estimated by summing the storage term, the resolved E C flux term (p w'c'), and the a  vertical SGS flux (Here after, we call it the 'mimic' E C approach).  Fig. 4.19 shows the comparison of the estimated N E E from the 'mimic' E C approach and the incorporated N E E . The strong effects of daytime convection have been largely eliminated by the mimic E C approach. Compared to the direct estimation of N E E without extracting half-hour averages (shown in Fig. 4.18b), the estimated N E E by the 'mimic' E C approach get fairly close to the incorporated N E E for the full 24-hour period, and particularly at night the underestimation of nocturnal respiration is within 20% of the prescribed source value, which reflects the contribution of horizontal advection at night. These two figures (Fig. 4.18b and Fig. 4.19) show very clearly the utility of the model for helping to establish and/or validate data-analysis procedures. However, although the estimated N E E by the 'mimic' E C approach is within 20% of the real source the vertical fluxes are dominated by the sub-grid-scale vertical flux. Thus, further evaluation is needed using much higher model resolution to see i f this result still holds when vertical fluxes are dominated by resolved-scale fluxes. Also, the dramatic positive effect of subtracting temporal averages suggests that further numerical studies could be used to assess optimal averaging techniques to further improve the E C method. It seems quite likely that such results would be site dependent. In summary, the contribution of each term in the CO2 budget varies with time (Fig. 4.13). Taking only horizontal or vertical advection into account to estimate the source, one would end up with an incorrect conclusion (Fig. 4.15). The integrations of each term in the CO2 budget from the ground to the height of the tower instruments (42 m) comprise the estimation of N E E . The storage term (Fig. 4.16) has near-zero value during daytime and rapid accumulation around sunset and quick exhalation around sunrise. The horizontal and vertical advection have opposite signs at day and largely offset each other, but they have same positive signs at night, indicating C 0 2 is transporting away from the tower site both by horizontal and vertical advection. At the West Coast Flux station, the sloping terrain contributes to the potential for thermally driven quasi-stationary  mesoscale  circulations (such as up- and down-slope flows, and land/sea breezes), resulting in significant horizontal and vertical advection. Under fair weather, the mean vertical velocity associated with convective turnover could also result in significant vertical advection.  horizontal advection (term [2]+[4J) vertical advection (term [3]+[5])  30  Time ( h i # 12 horizontal advection (term [2]+[4]) vertical advection (term [3]+[5])  10 eg  e  8  6 O E 4 o o 2 E 0 c o o -2 CD > -4 -o < -6 -8 20  21  22  23  Tfffie ( r ^ l )  26  27  28  29  Figure 4.17: Time series of the contributions of horizontal and vertical advection in NEE estimation (4-12). The lower panel shows a part of the upper panel from 20 to 29 h (corresponding to 8 pm to 5 am LST). The light line represents horizontal advection, and the dark line represents vertical advection.  40  estimated N E E b y the s u m of 5 terms incorporated real source (NEE)  30  ST  2o  CM  o  i  o  o 111 IXI-20 -30 -40  15_.  20  Time  -50  30  (hrs)  e s t i m a t e d N E E u s i n g t e r m [1-+-3+5] incorporated real s o u r c e  15  20  Time  25  30  (hrs)  Figure 4.18: Comparison of the incorporated source (thick line) to the estimated source using: a) the sum of five terms (light line in upper panel), and b) the sum of three terms (light line in lower panel). 40 30  Jj o  10  I  0  o  -H-10  UJ  2-20 -30 estimated NEE using mimic E C approach incorporated real source  -40 -50  10  15  20  25  30  Time (hrs)  Figure 4.19: Comparison of the incorporated source (the thick line) to the estimated source using an analysis procedure that mimics the EC approach (the light line).  The best way to estimate N E E is to sum all five terms (Fig. 4.18a). Using the sum of storage term, resolved and SGS vertical fluxes to estimate N E E , without reducing the strong effects of convection, shows large fluctuation during daytime and underestimates nocturnal respiration by about 40% of prescribed N E E (Fig. 4.16; Fig. 4.18b). The estimated N E E by the 'mimic' E C approach (Fig. 4.19) gets fairly close to the incorporated N E E for the full 24-hour period, due to the reduction in the effects of daytime convection. At night, the underestimation of nocturnal respiration is within 20% of the correct value. However, the vertical fluxes are dominated by the sub-grid-scale term, which indicate that the resolution in the current simulations is still not high enough to solve the essential physics. Further evaluation is needed to see i f this result still holds when vertical fluxes are dominated by resolved-scale fluxes.  4.5 Summary and Conclusions This chapter describes the parameterization schemes we adopted for the treatment of the sources and sinks in the CO2 budget equation. They are incorporated into a mesoscale atmospheric model (WFIS) and solved in the context of a diurnal cycle over a two dimensional mountain of similar horizontal and vertical scale as Vancouver Island. A l l meteorological variables and the CO2 mixing ratio are output every 15 seconds from the simulations and are used to perform the CO2 budget analysis and the N E E estimations. The contribution of each term in the CO2 budget is investigated and the effects of advection on N E E estimates at the West Coast Flux station are evaluated. The simulated vertical profiles and time series of CO2 concentration show a significant diurnal variation near the surface, corresponding to the diurnal change of stability in the A B L . They agree well with some observations in the literature [Yang, et al., 1999; Y i , et al., 2000; Eugster and Sigrist, 2000]. CO2 is considerably well mixed during the day and rapidly accumulates in the trunk space around sunset, resulting in a sharp vertical C O 2 gradient. This CC^-rich layer deepens and reaches its maximum depth near midnight. For the remaining night the CO2 concentration slowly decreases over the land. In contrast, the CO2 profiles over the neighboring water surface show a continuous increase of C O 2 throughout the night in spite of no locally respired CO2 efflux. Visualizations of the time evolution of the C O 2 spatial field show that the locally respired CO2 on the land is  transported from the land to the water surface at night by mesoscale advection processes such as slope flow and the land breeze. This transport is consistent with the simulations of wind at the West Coast Flux station in Chapter 3. The simulated vertical CO2 gradients above 9 m are small but below 9 m they are large. These gradients are greatest around sunset and sunrise, when the accumulation of CO2 near the surface and the continuous loss of CO2 in the upper canopy result in a rapid increase of vertical gradients. The time of occurrence of the sharp vertical CO2 gradients is in good agreement with the occurrence of the anomalies shown in the flux measurements around sunrise and sunset (Fig. 1.3 and Yang, et al., 1999). The simulated horizontal CO2 gradients are much smaller than the vertical gradients. The horizontal gradients result mainly from varying soil temperature due to variations in elevation. The CO2 budget (4-11) analysis further quantifies the contributions of each term in the CO2 budget at the West Coast Flux station, and shows the significance of advection in assessing whether an ecosystem is a source or sink. The contribution of each term in CO2 budget varies with time. During daytime, the horizontal advection and vertical advection show significant fluctuation, have opposite signs and largely counteract each other (Figs. 4.13, 4.16, 4.17). The net effect is a near zero storage of CO2, responding to strong daytime eddy mixing. A t night, horizontal and vertical advection are both responsible for transporting  away  the  locally  respired  CO2  in  the  trunk  space  (w<Oandc3c7az<0, « > 0 a n d ^ / a y > 0 , see Figs. 3.10, 4.8, 4.9, 4.11, 4.17). Around sunset CO2 starts to accumulate in the trunk space and still loses in the upper canopy. The magnitude of horizontal SGS flux is the smallest one and it is reasonable to ignore this term. A l l other terms significantly contribute to the CO2 budget at the West Coast flux station. The simulations and CO2 budget analysis support the hypothesis that advection plays an important role in the CO2 budget at the West Coast Flux station. Thermally driven slope flow, convection and land/sea breezes are the major mechanisms advecting locally respired CO2 away. The integration of each term in the CO2 budget from the ground to the height of the tower instruments (42 m) comprises the estimation of the N E E budget (4-12). The integration without extracting temporal (half-hour) means from each field variable shows that the advection terms (both horizontal and vertical) are the dominant terms both during  the day and at night. Ignoring advection horizontally or vertically could result in underestimating N E E at night and losing sight of the real transport of CO2 during the day. When half-hour means are subtracted from each field (the 'mimic' E C approach) to estimate the N E E , consistent with the E C approach currently used at F L U X N E T sites, the estimated N E E is much closer to the prescribed source due to the elimination of strong convection during daytime. A t night, the underestimation of nocturnal respiration is also reduced, and is within 20% of the prescribed source value. However, the estimated N E E at night is dominated by the sub-grid-scale vertical flux. Thus, further evaluation is needed using much higher model resolution to see i f this result holds when vertical fluxes are dominated by resolved scale fluxes.  Exploratory analysis of the advection source area  5.1 Introduction To investigate the source region of air parcels that are advected in the A B L to the instrumented tower, meteorological streaklines and trajectories are useful [Apadula et al., 2003; Aalto et al., 2003; Stohl et al., 2002]. Streaklines are lines of a real or hypothetical tracer deposited in the flow during a time interval by continuous emission from a fixed point [Stull, 2000]. A streakline shows the resulting plume at one instant in time, such as a satellite photo. A backward streakline indicates the locus of air parcels that will eventually pass through a receptor point at some later time. Namely, it can be used to determine where these air parcels come from, and where they move en-route to the receptor point. Trajectories, also called path lines, trace the route traveled by an air parcel during a time interval [Stull, 2000]. Essentially, it is the time integration of the position of an air parcel as it is transported by the wind [Stohl et al., 2002]. The trajectories can be integrated both forward and backward in time to produce forward and backward trajectories. The destinations of CO2 molecules emitted from a known source can be found using a forward streak or trajectory along the mean wind, while the source locations of CO2 molecules arriving at a known receptor or tower site can be found from a backward streakline or a backward trajectory. In this study, backward streaklines for the two tower sites at the West Coast Flux station were calculated daily from the surface wind fields generated by a mesoscale real-time weather prediction model M C 2 , which is run operationally by the Weather Forecast Research Team (WFRT) at the University of British Columbia. The hourly data files generating these backward streaklines were collected daily and plotted to identify the mesoscale advection source areas. This is the first step toward eventually building a sourcefootprint climatology for the tower sites. Due to the time limitations, the exploratory analysis was performed for only one season from the late fall to early winter (Oct., Nov. and Dec.) in 2004. Thus, this part of the study shifts away from the idealized fair-weather assumptions of the previous chapters, and instead considers the normal parade of stormy  and fair-weather synoptic and mesoscale events that make up typical west coast weather in British Columbia.  5.2 Methodology Three components comprise this source area modeling strategy: 1) Mesoscale modeling of the surface wind fields from the M C 2 model; 2) Application of the P U F F model to the surface wind fields generating the backward streaklines and trajectories; and 3) Daily data collecting and analyzing of the backward streaklines to locate the source positions of the air parcels that will pass through the flux receptor at two specified times. The first step is the daily generation of surface wind fields covering the flux towers, which was done by Xingxiu Deng, who was responsible for operationally running the weather prediction model — M C 2 . The M C 2 model is running with five-nested domains, where grid resolutions are 108 km, 36 km, 12 km, 4 km and 2 km, respectively (Fig. 5.1). The 2 km domain encompasses the Mature Stand tower site near Campbell River (Northern site, 49.905°N, 125.366°W) and the Young Fir tower site (Southern site, 49.52°N, 124.9°W). The surface wind fields are generated daily from the M C 2 2-km-grid real-time run, which is initialized at 1200 U T C and runs for a 27 hour forecast period.  125W Figure 5.1: MC2 model domain configuration. Left panel shows the 108 km, 36 km, 12 km and 4 km domains. The right panel shows the 4 km and 2 km domains. (Provided by Xingxiu Deng) The second step is the daily generation of the back streaklines and back trajectories  using the surface wind for the two flux tower sites, which is done using a program called P U F F written by George Hicks II in the W F R T . In order to construct the location and the pathways of air parcels arriving at the CO2 monitoring sites for the whole forecast period, the backward streaklines and trajectories for each hour are computed by the P U F F model using the surface wind fields (10 m above ground) provided by the M C 2 model, so these are not 3D trajectories. The basic idea of the P U F F model is that the program takes near surface winds and does a spatial and temporal linear interpolation along each time-step to determine the location and the movement of the air parcels. Dispersion and vertical movements are not considered in the P U F F model. We track back the positions of every air parcel that will eventually pass through the tower site and plot their locations once every hour on a separate "map". For a 27-hour forecast period of the 2-km-grid run, there are 27 maps showing the locations and trails of the air parcels in each hour. A l l the maps can be combined into an animation (movie) that runs for about 27 hours into the future, which are shown daily on the website: [http://weather.eos.ubc.ca/wxfcst/]. A n example map (Fig. 5.2) is arbitrarily chosen showing back streaklines with back trajectories and surface wind (10 m) at the Mature Stand flux site from our database. The back streaklines in the figure look like thick lines that actually consist of many colored dots, representing the air parcels that pass through the flux tower every 10 minutes. The dot colors change every 2 hours. B y associating the time with the corresponding color of the parcels (indicating by the color bar), we can easily track the locations and the movement of different parcels in a backward streakline, thus to locate where an air parcel comes from, and where it moves en-route to the tower. The thin lines in Fig. 5.2 represent the backward trajectories that are the trails left by the air parcels in the back streakline. The backward trajectory of an air parcel that passed over the tower at one time is usually different than the backward trajectory for a different time. Hence, the air parcels leave different trails (i.e., different back trajectories). The net result is an advective "footprint", showing all regions of the surface traversed by each air parcel as it travels toward the tower. The final step is the daily data collection and analysis of back streaklines to locate the source positions of the air parcels passing over the flux tower, which was my task. The flux tower researchers are mainly interested in the source position of air parcels 8 hours before  Model: UC2  initioiized: OUTC, 29 DEC 2 0 0 4  Surface Wind Direction and Back Streaks  2 km grid  Forecast valid 16PST 29 DEC (17MST 29 DEC) [OOUTC 3D DEC] 2004  0  2  4  6  S  10  12  14  16  18  20  22  24  colour bar: Time an air parcel reaches the tower in UTC arrows: Wind Direction  Figure 5.2: Surface wind (arrows) and back streaks (thick line consisted of colored dots) with trajectories (thin lines) passing through the northern flux site at 16 pm, Dec. 29, 2004. The time of day (in PST, MST, and UTC) is shown at both the top and bottom of the map. The color bar indicates the UTC time when the same color parcel pass through the tower site.  passing through the tower at 4 pm and 4 am local standard time (LST) (corresponding to 2400 and 1200 UTC). Every day, we located the exact position of the air parcel that will pass through the tower site at those two times. The source locations of these air parcels have been collected daily for three months: October, November and December 2004, for both the Mature stand site (northern site) and Young Fir site (southern site) at the West Coast Flux station. The results are shown in Fig. 5.3, 5.4, 5.5 and 5.6.  5.3 Results Figure 5.3 shows the source locations of the air parcels 8 hours before they passed through the Mature Stand site at 1600 L S T in October (a), November (b) and December (c), and Oct. to Dec. (d) in 2004. In Oct. and Nov. for those parcels arriving the tower at 1600 L S T (here after, we call them 'daytime arriving parcels'), their sources were mainly located along Georgia Strait, but there were two incoming channels: from the northwest and from the southeast, corresponding to the geographic orientation of Georgia Strait. The parcels coming from the southeast generally traveled a longer distance along Georgia Strait to eventually reach the tower perhaps because of less drag on the water surface allowing faster winds. Thus, the source positions of these parcels were more sparsely distributed along Georgia Strait. In D e c , the source positions for the daytime arriving parcels are densely gathered around the tower because those parcels mainly came from the northwest, where the wind was generally slower. The combined plot (Fig. 5.3 (d)) further shows for the whole season (from Oct. to Dec.) that the source positions of daytime arriving parcels at the Mature Stand site were mainly distributed from the southeast to northwest along Georgia Strait, covering an area of thousands of km . The source area of CO2 concentration at the 2  Mature Stand site was greatly influenced by the northwesterly and southeasterly flow along Georgia Strait. Fig. 5.4 is same as Fig. 5.3 except for those parcels passing through the Mature Stand site at 0400 L S T (here after, we call them 'nighttime arriving parcels'). During nighttime, the wind speeds were generally lower than at day and the paths of the parcels traveling were much shorter during the 8-hour period. The source locations of the nighttime arriving parcels were densely distributed around the tower site, except 9% outliers spreading over Georgia Strait. For daytime arriving parcels, more than 50% of them had source locations  over Georgia Strait, even extending to the Greater Vancouver area (Fig. 5.3). Therefore, compared to daytime arriving parcels, the source area of the nighttime arriving parcels was much smaller. That suggests that during nighttime, with typically low wind speed, the CO2 concentration in the shallow stably stratified A B L is dominated by the local CO2 flux, which is controlled by the local downslope drainage flow discussed in earlier chapters. However, during the daytime with deep well-mixed A B L , the CO2 concentration could be strongly influenced by mesoscale and synoptic-scale atmospheric transport. Figure 5.5 and 5.6 shows the source locations of air parcels 8 hours before passing through the Young Fir site at 1600 L S T and 0400 L S T , respectively. The Young Fir site is much closer to Georgia Strait compared to the Mature Stand site, and thus is affected more by the prevailing wind over Georgia Strait. Fig.5.5 shows that more than 70% of the daytime arriving parcels came from Georgia Strait over the water, mainland and islands and Vancouver — very heterogeneous surfaces. The nighttime arriving parcels are similar to those from the Mature Stand site, namely their source locations are generally not too far from the Young Fir destination. Thus, the advection influence on the CO2 concentration measurements mainly came from the local region. Fig. 5.5 (d) and 5.6 (d) show that the source areas of the Young Fir site during daytime (1600 LST) were mainly along Georgia Strait and for nighttime (0400 LST) are mostly local for the whole season.  5.4 Conclusion The CO2 concentration measurements generally have much larger footprints than eddy flux measurements because of the influence of advection, which is also shown in footprint analysis [e.g. Kljun et al., 2002; Schmid, 1994]. To aid the interpretation of the C O 2 concentration measurements observed at the tower sites and to better understand the scaling of the local CO2 measurements to regional scales, back streaklines at two flux tower sites were used to identify the advection source area. The results show that for Oct., Nov., and Dec. 2004, the daytime arriving parcels mostly came from the south portion of Georgia Strait and traveled along Georgia Strait for over 200 km during the 8-hour period. The advection source area during daytime at the West Coast Flux station mainly extended along Georgia Strait and covered over thousands of km , corresponding to the prevailing 2  northwesterly and southeasterly flows over Georgia Strait. The nighttime arriving parcels  had a very limited source area, coming from the surroundings of the tower sites. Therefore, we conclude that at nighttime with typically low wind speed, the CO2 concentration in the shallow stably stratified A B L is dominated by the local CO2 flux. During the daytime with deep well-mixed A B L , the C O 2 concentration is strongly influenced by mesoscale atmospheric transport. Comparing the results of this chapter with the 2D simulations from the earlier chapters, it appears that the 2D simulations across Vancouver Island do not adequately represent the near-surface winds that are strongly channeled by Georgia Strait. However, a related criticism is that the back-streak model utilized only the surface winds and not the full 3-D wind field. A n obvious recommendation for future work is that both the WFIS large-eddy simulations and the P U F F real time trajectory forecasts must be expanded to their full 3-D capabilities to adequately capture the complex flow near the flux tower sites. ~  50.4N  ~ 0  -\.' : -  -  .'•"''/  MoUire Stand  1  ooe ' • 00  SON 49.8N  49.8N  49.6N 49.4N  Mclure Stcnd Nov. 16pni  0  50.2N  0  1°  0 0 v .  49.6N  \  °° . ° V.  / i  %o  0  124W  123W  0 : :  0  o 0  122W  49N 126W  125W  ':  0  49.2N  a 125W  -.  49.4N  1  49.2N 49N !26W  0 0  0  I24W  b  123W  122W  Mature Stand Oct.-Dsc. 16pm  0 O  I24W  123W  122W  49N-! 126W  125W  124W  0  ; - 0  123W  122W  Figure 5.3: The locations of air parcels (circles) 8 hours before passing through the Northern site (triangle) at 1600 LST in October (a), November (b) and December (c), and Oct. to Dec. (d), 2004.  50:4N  50.4N  Molure Stend Oct. 4am  50.2N 50N 49.8N 49.6N  50N •J * o  oo  \ tf>  49.8N ••. ' •••  .,  49.6N  0  49.4N  ....  v  s  49.4N  49.2N  49.2N  a  49N 126W  125W  124W  50.4N  123W  122W  49N 126W  I24W  SON  50M 49.8N  123W  122W  Maiure Stcnd Oct.-Dec. 4am  50.2N  49.8N °  b 12SW  50.4N  Maxure Stond Dec. 4am  50.2N  o  49.6N  49.6N  49.4N  49.4N  49.2N  49.2N  49N-I 126W  Molure Stond Nov. 4am  50.2N  125W  124W  123W  122W  49N 126*  125W  124W  123W  122W  Figure 5.4: The locations of air parcels (circles) 8 hours before passing through the Northern site (triangle) at 4 am PST in October (a), November (b) and December (c), and Oct. to Dec. (d) 2004. 50.5N  50.5N Younq Fir 16pm  QrA.  50H  Younq Fir  iiov. 16pm  50N  49.5N  49.5N  49N  49N  o 48.5N  48N-1 126W  .-o  48.5M 125W  124W  123W  122W  50.5N  48N 126W  125W  I24W  123W  122W  50.5N o ;,>,.•.., ! 50N-  49.5N  Younq Fir Oct.-0ec. 16pm  ° % o aO c °o  :\o 9 *'<>,,,  49N  48.5N  48N 126W  125W  124W  123W  122W  4BN 126W  125W  124W  I23W  Figure 5.5: Same as Fig.5.3 except for the Southern (Young Fir) site.  122W  50.5N ,  Younq Fir Oci. 4om  0  <&> 0  ©  0  0 0  a  0  0 0  a 126W  125W  124W  123W  122W  126W  125W  124W  123W  122W  50.5N s  Younq Fir Dec. 4 am  Younq Fir Oct.-Dec. 4am  50NH 0  ft / 0  49.5NH  0 "  00  o «  0°  ..  49NH • - \  (O 0  ,' -  <fi  • 0  0 0  0 0  48.5N H  c 126W  125W  124W  I23W  122W  48N ~i'26W  125W  124W  123W  Figure 5.6: Same as Fig. 5.4 except for Southern (Young Fir) site.  I22W  Summary, Conclusions, Limitations and Future work 6.1 Summary and Conclusions The E C approach used to estimate N E E (storage plus eddy flux) ideally requires that a flux tower be located on flat and homogenous terrain so that the mean horizontal and vertical advection effects are neglected, and the effects of the mean horizontal turbulent flux are small compared to the vertical. However, sloping and heterogeneous topography is found instead at some flux tower sites, which could strongly affect the application of E C method on estimating the strength of carbon sequestration or loss. The West Coast Flux station has such complex topography, being located on the lee side of Vancouver Island on a 5-10 degree slope. Georgia Strait is only 9 km away from the tower site. The complex (coastal and mountainous) topography raises uncertainties in the E C estimation of N E E because of the underlying advective influence of mesoscale circulations such as slope flows, land/sea breezes and convection. However, these mesoscale mechanisms are much more difficult to measure. Motivated by the reported uncertainties and difficulties in observation, this study employed a high-resolution atmospheric model, WFIS, to investigate the roles of mesoscale flow on the E C measurements, and to evaluate the advection effects on N E E by quantifying the contributions of each term in the CO2 budget. In particular, the goal of the research is to aid in the understanding of the mesoscale controlling mechanisms that might be responsible for transporting C O 2 away from the flux tower site. This could help account for the anomalies in eddy-covariance flux measurements at the West Coast Flux station. To enable this study, WFIS was refined to treat canopy flow and CO2 fluxes. The refinements of WFIS included (a) adding the tree-drag force in the momentum equations; (b) incorporating long-wave and short-wave radiation effects of the canopy on the surface and canopy energy budget; (c) adding a CO2 budget into the model with a treatment of source/sink terms such as canopy photosynthesis and soil respiration; and (d) incorporating a four-layer soil model using a heat transfer equation to simulate the soil temperatures at different depths.  Results and conclusions are: 1.  The experiments of sensitivity to drag coefficient (Appendix C) suggest that an deep forest canopy increases the tendency for the lee flow to separate and reverse i f the tree-drag is sufficiently large, i.e., a sufficiently dense canopy. The tests of sensitivity to upstream flow and topography indicate that higher hills result in stronger flow impinging on the canopy, and lead to earlier flow separation and stronger flow reversals at the downstream edge of the canopy.  2.  A dense canopy reduces the diurnal range of surface air temperature and soil temperature. The surface air and soil temperatures under a dense canopy are uniformly lower than over a bare soil surface during the day and higher at night (Figs. 3.3, 3.4, 3.11).  3.  The simulated winds at the West Coast Flux station show a distinct diurnal pattern above the canopy. During the day, the winds flow upslope due to the anabatic flow and sea breezes, at night the flows are dominated by down-slope flow and weak land breezes (Figs. 3.5, 3.6, 3.8, 3.9). But under the canopy, the statistics of the simulated winds indicate that 60% of the time during a day downslope flows occur because of the convective thermal turnover (Fig. 3.10). The wind pattern around sunset shows a quick transition of wind direction and much weaker wind speed (Fig. 3.6). The rather rapid shift from upslope to down-slope winds just after sunset may be an important factor in helping explain some of the CO2 flux anomalies just after sunset at the West Coast Flux tower.  4.  The drainage flow for the dense-canopy case is much weaker and shallower compared to the bare-ground case, and is separated into sub-canopy and abovecanopy parts by the temperature inversion at the canopy top (Fig. 3.7). The wind profiles indicate that the minimum wind occurs in the upper canopy, and then increases with decreasing height above the ground, until surface friction reverses this trend close to the ground level (Figs. 3.12, 3.13).  5.  The simulated vertical profiles (Fig. 4.6) and time evolution of CO2 concentration (Figs. 4.8, 4.9, 4.10) show a significant variation of CO2 near the surface, corresponding to the diurnal change of stability in the A B L . The CO2 field is considerably well mixed during the day, and rapidly accumulates in the trunk  space around sunset. The COvrich layer deepens and reaches its maximum depth near midnight. For the rest of the night the CO2 concentration slightly decreases over the land. 6.  The simulated vertical CO2 gradients above 9 m are small, but below 9 m they are large (Fig. 4.8, 4.9). This is particularly the case around sunset and sunrise, when the accumulation of CO2 near the surface and the continuous loss of CO2 in the upper canopy cause a rapid increase of vertical gradients. The time of the sharp CO2 gradient occurring is in good agreement with the time of CO2 flux anomalies shown in the measurements (Fig. 1.3 and Yang et al., 1999). The simulated horizontal CO2 gradients are much smaller than the vertical gradients (Fig. 4.11). The horizontal gradients result mainly from the varying soil temperature due to variations in elevation.  7.  CO2 advection is indicated by the accumulation of CO2 concentration over the neighboring water surface at night (Fig. 4.7). Visualizations of the time evolution of the CO2 spatial field (Fig. 4.12) show that the locally respired CO2 on the land is transported to the neighboring water (Georgia Strait) at night by mesoscale mechanisms such as the down-slope flows and land breezes under fair-weather conditions at the West Coast Flux station.  8.  The CO2 budget analysis further quantifies the contributions of each terms in the CO2 conservation equation, and supports the hypothesis that advection plays an important role in the CO2 budget at the West Coast Flux station. The contribution of each term in the C O 2 budget varies with time. During daytime, vertical advection is generally responsible for bringing CO2 into the control volume, while horizontal advection usually accounts for CO2 loss (Figs. 4.13, 4.14). The net effect is near-zero storage of CO2 within the canopy space responding to the strong eddy mixing. At night, horizontal advection and vertical advection are both responsible for transporting the locally respired CO2 away near the surface (Figs. 4.13, 4.14). The magnitude of the horizontal SGS flux is the smallest, and it is reasonable to ignore this term. A l l other terms are significant contributors to the simulated CO2 budget at the West Coast flux station.  9.  The integration of each term in the CO2 budget from the ground to the height of the tower instruments (42 m) comprises the contribution of each term to the estimation of the N E E budget (4-12). The storage term (term [1] in Fig. 4.16) has near-zero value during daytime and rapid accumulation around sunset and quick exhalation around sunrise. The horizontal and vertical advection have opposite signs at day and largely offset each other, but they have same positive signs at night, indicating C 0 2 is transporting away from the tower site both by horizontal and vertical advection (Fig. 4.16, 4.17). Using only the sum of the storage and vertical-flux terms to estimate N E E , without reducing the effects of convection, resulted in strong fluctuations during the day and the underestimation of nocturnal respiration by about 40% (Fig. 4.18b). When we subtracted half-hour averages (as is done in the E C technique), the daytime fluctuations of estimated N E E are strongly reduced and the nocturnal respiration is within 20% of the prescribed source, indicating a 20% contribution to N E E by the horizontal advection (Fig. 4.19).  10. The back-streakline analysis suggests that in Oct., Nov., and Dec. 2004, the advection source area of CO2 concentration during daytime mainly extended along Georgia Strait and covered over thousands of k m northwest and southeast 2  of the Flux station (Fig. 5.3, 5.5). However, at night there was a smaller CO2 advective footprint closer to the West Coast Flux station (Fig. 5.4, 5.6). Therefore, we conclude that at nighttime with typically low wind speed, CO2 concentration in the shallow stably stratified A B L is dominated by the local CO2 flux and local advection, as was simulated by WFIS. During the daytime, with deep well-mixed A B L , the C O 2 concentration is strongly influenced by mesoscale atmospheric transport  from locations far outside of the WFIS  simulation domain. In short, the calculations examined the impacts of a deep forest canopy on the surface energy budget and the A B L flow, and assessed the advection effects on the N E E estimation. Over an idealized two-dimensional mountain of similar horizontal and vertical scale as Vancouver Island, the 2D vertical slice experiments capture the firstorder effects of diurnal heating/cooling on the sloping terrain. The volume exchange of  heat, momentum and CO2 above and within the canopy appears to be strongly affected by local flows resulting from diurnal thermal forcing. If synoptic and large mesoscale flows are neglected, then the anabatic/katabatic flow, land/sea breezes and convective thermals under fair weather conditions are the major mechanisms for advection of CO2 at the West Coast Flux station. The simulated C O 2 concentration shows significant variation and gradients near the surface, and captures the transport processes at night. The C O 2 budget analysis further quantifies the contributions of advection (both horizontal and vertical). The estimation of N E E , without reducing the effects of strong convection, shows large fluctuations during daytime, and underestimates  nocturnal  respiration by about 40%. With the subtraction of the half-hour mean, the estimated N E E is much closer to the prescribed source. Back-streakline analysis provides insight on the advection source areas of CO2 concentration when synoptic and mesoscale flows typical of fall weather (both fair and stormy) are not neglected. Substantial advection parallel to Georgia Strait is anticipated during the other seasons too, but was not examined yet.  6.2 Limitations A l l numerical experiments and analysis in the current research were performed using twodimensional (2D) vertical across-section applications of the WFIS model, and the quasi-2D horizontal (surface-following) winds of M C 2 model for the back-streak analysis. These 2D simulations and analysis provided theoretical and physical insights for understanding the real three dimension (3D) atmospheric conditions, and were also useful for testing the new parameterization schemes we incorporated into the WFIS model. However, some of the limitations of 2D simulations are obvious. The first drawback of the 2D experiments is that we assume that there is no lateral motion that can transport CO2. The horizontal advection effects we investigated with WFIS were only along the slope. Thus, the 2D experiments cannot capture ridge-parallel advection effects that appeared to be prominent in the back-streak analysis using M C 2 surface wind. The second shortcoming of the 2D WFIS experiments is that the dynamics in these experiments has only 2D eddy characteristics (upscale nonlinear energy cascade) rather than 3D characteristics (downscale nonlinear energy cascade). This effect would probably be strongest during the daytime periods of highly nonlinear convection. At this point, the robustness of some explanations  in this study such as the explanation of the rapid transition of the wind pattern around sunset will still need to be tested using 3D simulations. The vertical SGS term dominates at night (Fig. 4.16, 4.19), indicating the resolution we used in this study (75 m in horizontal and 6 m in vertical) is still not high enough to resolve the essential physics within the canopy. Another limitation is the assumption of zero stemdrag in the trunk space, which may allow the wind to be too fast within the trunk space, causing excessive transport of CO2. The assumption of no heat storage within the canopy may also somewhat affect the simulations of the surface energy budget.  6.3 Future work The current research could be extended to include the following: 1.  The limitations of the 2D WFIS framework could be resolved by considering a 2D3D experiment [e.g. Clark and Farley, 1984], where a limited number of grid points are added in the third dimension to allow some 3D dynamics.  2.  The current 2D parameterization schemes that we added in the model could be modified to work in 3D. The 3D application of WFIS could be utilized to test the robustness of those results gained from the 2D experiments in the current research.  3.  The current experiments used an idealized two-dimensional cross-section hill representative of Vancouver Island through the tower site for the topography. In future 3D experiments, we could take the full 3D topography into consideration within a full 3D WFIS model.  4.  In the current simulations, we used a homogenous canopy distribution on the land to simulate the influence of a deep forest canopy on the A B L flow. In the future, we could use a more realistic canopy distribution that includes the effects of varying canopy height in the future.  5.  Some real-time cases could be run to compare the simulated dynamic and thermodynamic fields with real observations. Analysis of the CO2 budget could be performed using output from the full 3D experiments to investigate advection effects on the flux measurements at the West Coast Flux station.  6.  So far, the back-streaklines analysis in the current study is only based on the data for three months. The back-streak analysis could be extended to include the data in one or two years of data to build a full climatology of the advection source area for the tower sites at the West Coast Flux station. Also, the P U F F program could be modified to compute the full 3D back streaks.  7.  We could further improve radiation budget to adopt other solar radiation option in WFIS model to include the diffuse short wave radiation under cloudy condition.  8.  We could treat tree drag as a function of Froude number for different atmospheric stability.  9.  In this research the vertical sub-grid-scale CO2 flux dominates at night, which indicate that the resolution we used is still not high enough to solve the essential physics, we could consider higher resolution in the future.  10.  When we perform the C 0 2 budget analysis by mimicking the E C approach, the density correction in Web et al. [1980] could be considered in the future.  Bibliography Aalto, T., J. Hatakka and Y. Viisanen, 2003: Influence of air mass source sector on variations in C 0 mixing ratio at a boreal site in northern Finland. Boreal Envir. Res., 8, 385-393. 2  Allen, T. and A.R. Brown, 2002: Large-eddy simulation of turbulent separated flow over rough hills. Boundary-Layer Meteorol., 102, 177-198. Amiro, B.D., 1990: Drag coefficients and turbulence spectra within three boreal forest canopies. Boundary-Layer Meteorol., 52, 227-246. Anthoni, P.M., B.E. Law and M.H. Unsworth, 1999: Carbon and water vapor exchange of an opencanopied ponderosa pine ecosystem. Agric. and Forest Meteorol, 95, 151-168. Apadula, F., A. Gotti, A. Pigini, A. Longhetto, F. Rocchetti, C. Cassardo, S. Ferrarese and R. Forza, 2003: Localization of source and sink regions of carbon dioxide through the method of the synoptic air trajectory statistics. Atmo. Envir., 37, 3757-3770. Arakawa, A., 1966: Computational design for long term integration of the equations of fluid motion: two-dimensional incompressible flow. Part I. J. Comput. Phys., 1, 119-143. Arya, S. P., 1988: Introduction to micrometeorology. Academic Press, Inc. 307 pp. Atwater, M.A. and P. S. Brown, 1974: Numerical calculation of the latitudinal variation of solar radiation for an atmosphere of varying opacity. J. Appl. Meteorol, 13, 289-297. Ayotte, K.W., J.J. Finnigan and M.R. Raupach, 1999: A second-order closure for neutrally stratified vegetative canopy flows. Boundary-layer Meteorol, 90, 189-216. Aubinet, M . , A. Grella, A. Ibrom, U . Rannik, J. Moncrieff, T. Foken, A. Kowalski, P. Martin, P. Berbigier, C. Bernhofer, R. Clement, J. Elbers, A. Granier, T. Griinwald, K, Morgenstern, K . Pilegaard, C. Rebmann, W. Snijders, R. Valentini and T. Vesala, 2000: Estimates of the annual net carbon and water exchange of European forests: the Euroflux methodology. Advances in Ecological Research, 30: 113-175. Baldocchi, D.D., B.B. Hicks and T.P. Meyers, 1988: Measuring biosphere-atmosphere exchange of biologically related gases with micrometeorological methods. Ecology, 69(5), 1331-1340. Baldocchi, D. and T. Meyers, 1998: On using eco-physiological, micrometeorological and biogeochemical theory to evaluate carbon dioxide, water vapor and trace gas fluxes over vegetation: a perspective. Agric. For. Meteorol, 90: 1-25. Baldocchi, D.D., E. Falge, L. Gu, R. Olson, D. Hollinger, S. Running, P. Anthoni, C. Bernhofer, K. Davis, J. Fuentes, A. Goldstein, G. Katul, B. Law, X. Lee, Y. Malhi, T. Meyers, J.W. Munger, W. Oechel, K . Pilegaard, H.P. Schmid, R. Valentini, S. Verma, T. Vesala, K. Wilson and S. Wofsy, 2001: FLUXNET: A new tool to study the temporal and spatial variability of ecosystem-scale carbon dioxide, water vapor and energy flux densities. Bull. Amer. Meteorol. Soc, 82, 2415-2435. Baldocchi, D.D., 2003: Assessing the eddy covariance technique for evaluating carbon dioxide exchange rates of ecosystems: past, present and future. Global Change Biology, 9, 479-492.  Belcher, S.E. and J.C.R. Hunt, 1998: Turbulent flow over hills and waves. Ann. Rev. Fluid Mech., 30, 507-538. Black, T.A., G. den Hartog, H.H. Neumann, P.D. Blanken, P.C. Yang, C. Russell, Z. Nesic, X . Lee, S.G. Chen, R. Staebler and M.D. Novak, 1996: Annual cycles of C 0 and water vapor fluxes above and within a Boreal aspen stand. Global Change Biology, 2, 219-229. 2  Bradley, F., 1980: A n experimental study of the profiles of wind speed, shearing stress and turbulence at the crest of a large hill. Q.J.R.M.S., 106, 101. Brown, A.R., J.M. Hobson and N . Wood, 2001: Large-eddy simulation of neutral turbulent flow over rough sinusoidal ridge. Boundary-Layer Meteorol, 98, 411 -444. Brunet, Y., J.J. Finnigan and M.R. Raupach, 1994: A wind tunnel study of air flow in waving wheat: single point velocity statistics. Boundary-layer Meteorol., 70, 95-132. Carrara, A., S.A. Kowalski, J. Neirynck, LA. Janssens, J.C.Yuste and R. Ceulemans, 2003: Net ecosystem C 0 exchange of mixed forest in Belgium over 5 years. Agric. For. Meteorol., 119, 209227. 2  Chen, W.J., T.A. Black, P.C. Yang, A.G. Barr, H.H. Neumann, Z. Nesuc, P.D. Blanken, M.D. Novak, J. Eley, R.J. Ketler and R.Cuencas, 1999: Effects pf climate variability on the annual carbon sequestration by a boreal aspen forest. Global Change Biology, 5 (1), 41-53. Cionco, R.M., 1965: A mathematical model for airflow in a vegetative canopy. J. Appl. Meteorol., 4, 517-522. Ciais, P., P.P. Tan and J.C. White, 1995: A large northern hemisphere terrestrial C 0 sink indicated by the C / C ratio of atmospheric C 0 . Science, 269, 1098-1102. 2  13  12  2  Clark, T.L., 1977: A small-scale dynamic model using a terrain-following coordinate transformation. J. Comp. Phys., 24, 186-215. Clark, T.L., 1979: Numerical simulations with a three-dimensional cloud model: Lateral boundary condition experiments and multicellular severe storm simulations. J. Atmos. Sci., 36, 2192-2215. Clark, T.L., 1982: Cloud modeling in three spatial dimensions. Hailstorms of the central high plains, Vol. 1, National hail research experiment, Charles A. K . and Patrick Squires, Eds. Colorado Associated University Press, 282 pp. Clark, T.L. and W.R. Farley, 1984: Severe downslope windstorm calculations in two and three spatial dimensions using anelastic interactive grid nesting: A possible mechanism for gustiness. J. Atmos. Sci., 41: 329-350. Clark, T.L. and W. D. Hall, 1991: Multi-domain simulations of the time dependent Navier Stokes equation: Benchmark Error analyses of nesting procedures. J. Comp. Phys., 92, 456-481. Clark, T.L. and W.D. Hall, 1996: On the design of smooth, conservative vertical grids for interactive grid nesting with stretching. J. Appl. Meteor., 35, 1040-1046.  Clark, T.L, W.D. Hall, R.M. Kerr, D. Middleton, L. Radke, F.M. Ralph, P.J. Neiman and D. Levinson, 2000: Origins of aircraft-damage clear-air turbulence during the 9 December 1992 Colorado downslope windstorm: Numerical simulations and comparison with observations. J. Atmos. ScL, 57, 1105-1131. Clark, T.L., M . Griffiths, M.J. Reeder and D. Latham, 2003: Numerical simulations of grassland fires in the Northern Territory, Australia: A new subgrid-scale fire parameterization. J. Geo. Res., 108, 1-15. Cox, P.E., A.R. Betts, A.D. Jones, S.A. Spall and I.J. Totterdell, 2000: Acceleration of global warming due to carbon-cycle feedbacks in a coupled climate model. Nature, 408, 184-187. Davidson, E.A., K. Savage, L.V. Verchot and R. Navarro, 2002: Minimizing artifiacts and biases in chamber-based measurements of soil respiration. Agric. For. Meteorol., 113, 21-37. Devonec, E. and A.P. Barros, 2002: Exploring the transferability of a land-surface hydrology model. J. Hydro., 265, 258-282. Drewitt, G.B., T.A. Black, Z. Nesic, E.R. Humphreys, E.M. Jork, R. Swanson, G.J. Ethier, T. Griffis and K. Morgenstern, 2002: Measuring forest floor CO2 fluxes in a Douglas-fir forest. Agric. For. Meteorol, 110 (4), 299-317. Eugster, W. and F. Siegrist, 2000: The influence of nocturnal C 0 advection on C 0 measurements. Basic Appl. Ecol., 1, 177-188. 2  2  flux  Fan, S., M . Gloor, J. Mahlman, S. Pacala, J. Sarmiento, T. Takahashi and P. Tans, 1998: A large terrestrial carbon sink in North America implied by atmospheric and oceanic carbon dioxide data and models. Science, 282, 442-446. Finnigan, J.J., 1979: Turbulence in waving wheat, II: Structure of momentum transfer. BoundaryLayer Meteorol., 16,213-236. Finnigan, J.J. and Y. Brunet, 1995: Turbulent air flows in forests on flat and hilly terrain. In Wind and Trees. Pp 3-40. Eds. M.P. Coutts and J. Grace. Cambridge University Press. UK. Finnigan, J.J., 2000: Turbulence in plant Canopies. Annu. Rev, Fluid Mech., 32, 519-571. Finnigan, J.J. and S.E. Belcher, 2002: Flow over a hill covered with a plant canopy. Q.J.R.M.S., 130, 1-29. Fitzjarrald, D.R., 1986: Slope winds in Veracruz. J. Clim. Appl. Met, 25, 133-144. Gao, W., R.H. Shaw, and K.T. Paw U, 1989: Observation of organized structure in turbulent flow within and above a forest canopy. Boundary-Layer Meteorol., 47, 349-377. Gal-Chen, T. and R.C.J. Sommerville, 1975: On the use of a coordinate transformation for the solution of the Navier-Stokes equations. J. Comput. Phys.,\l, 209-228. Garratt, J.R., 1992: The Atmospheric Boundary Layer. Cambridge University Press, 316 pp.  Grabowski, W.W. and T.L. Clark, 1991: Cloud-environment interface instability: Rising thermal calculations in two spatial dimensions. J. Amos. Sci. 48, 527-546. Grace, J., J. Lloyd, J. Mclntire, A.C. Miranda, P. Meir, H.S. Miranda, C. Nobre, J. Moncrieff, J. Massheder, Y . Malhi, I. Wright and J. Gash, 1995: Carbon dioxide uptake by an undisturbed tropical rain forest in southwest Amazonia, 1992 to 1993. Science, 270, 778-780. Holland E.A., 1999: North American Carbon sink. Science, 283, 1815a. Hoinka, K.P. and T.L. Clark, 1991: Pressure drag and momentum fluxes due to the Alps. I: Comparison between numerical simulations and observations. Q.J.R.M.S. 117, 495-525. Humphreys, E.R., T.A. Black, G.J. Ethier, G.B. Drewitt, D.L. Spittlehouse, E.-M. Jork, Z. Nesic and N.J. Livingston, 2003. Annual and seasonal variability of sensible and latent heat fluxes above a coastal Douglas-fir forest, British Columbia, Canada. Agric. Forest Meteorol. 115(1-2): 109-125. Hunt, J. C. R., K.J. Richards and P.W.M. Brighton, 1988a: Stably stratified shear flow over low hills. Q.J.R.M.S., 114, 859-886. Hunt, J. C. R., S. Leibovich and K J . Richards, 1988b: The turbulent shear flow over low hills. Q.J.R.M.S., 114, 1435-1470. Inoue, E., 1963: On the turbulent structure of air flow within crop canopies. J. Meteorol. Soc. Japan, 41,317-325. Jassal, R.S;, T.A. Black, G.B. Drewitt, D. Gaumont-Guay, M.D. Novak and Z. Nesic, 2004: A model of the production and transport of CO2 in soil: Predicting the C 0 efflux from a forest-floor and influences of soil water content and temperature. Agric. For. Meteorol., 124, 219-236. 2  Janssens, I.A., A.S. Kowalski, B. Longdoz and R. Ceulemans, 2001: Assessing forest soil C 0 efflux: an in situ comparison of four techniques. Tree Physiology, 20, 23-32.  2  Kaimal, J.C. and J.J. Finnigan, 1994: Atmospheric Boundary Layer Flows. Oxford University Press, 289 pp. Kessler, E., 1969: On the distribution and continuity of water substance in atmospheric circulations. Meteor. Monogr. 32, 84 pp. Kljun, N , M.W. Rotach and H.P. Schmid, 2002: A three-dimensional backward lagrangian footprint model for a wide range of boundary-layer stratifications. Boundary-layer meteorol. 103, 205-226. Koenig, L.R. and F.W. Murray, 1976: Ice-bearing cumulus cloud evolution: Numerical simulation and general comparison against observations. J. Appl. Meteorol, 15, 747-762. Komatsu, H , N . Yoshida, H. Takizawa, I. Kosaka, C. Tantasirin, and M . Suzuki, 2003: Seasonal trend in the occurrence of nocturnal drainage flow on a forested slope under a tropical monsoon climate. Boundary-Layer Meteorol., 106, 573-592. Law, B.E., D.D. Baldocchi and P.M. Anthoni, 1999: Below-canopy and soil C 0 fluxes in a ponderosa pine forest. Agric. For. Meteorol., 94, 171-188. 2  Lane, T.P., R.D. Sharman, T.L. Clark and H.M. Hsu, 2003: An investigation of turbulence generation mechanisms above deep convection. J. Atmos. Sci., 60, 1297-1321. Lavigne, M.B., M.G. Ryan, D.E. Anderson, D.D. Baldocchi, P.M. Crill, D.R. Fitzjarrald, M.L. Goulden, S.T. Gower, J.M. Massheder, J.M. McCaughey, M . Rayment and R.G. Striegl, 1997: Comparing nocturnal eddy covariance measurements to estimates of ecosystem respiration made by scaling chamber measurements at six coniferous boreal sites. J. Geophys. Res., 102 (28), 977-985. L i , Z.J., R. Miller and J.M. Lin, 1985: A first-order closure scheme to describe counter-gradient momentum transport in plant canopies. Boundary-Layer Meteor., 33, 77-83. Lee, x., T.A. Black, G. den Hartog, H.H. Neumann, Z. Nesic and J. Olejnik, 1996: Carbon dioxide exchange and nocturnal processes over a mixed deciduous forest. Agric. For. Meteorol., 81, 13-29. Lee, X . , 1998: On micrometeorological observations of surface-air exchange over tall vegetation, Agric. For. Meteorol., 91, 39-49. Lilly, D.K., 1962: On the numerical simulation of buoyant convection. Tellus, 14: 145- 172. Lipps, F.B., and R.S. Hemler, 1982: A scale analysis of deep moist convection and some related numerical calculations. J.Atmos. Sci., 39, 2192-2210. Lloyd, J. and J.A. Taylor, 1994: On temperature dependence of soil respiration. Func. Ecol, 8, 315323. Mahrt, L., 1982: Momentum balance of gravity flows. J. Atmos. Sci., 39, 2701-2711. Mahrt, L., 1986: Nocturnal topo-climatology. World Climate Applications Program Report WCP117,76 pp. Mahrt, L., 1998: Flux sampling errors for aircraft and towers. J. Atmos. Oceanic Techno!., 15, 416429. Manins, P.C. and B.L. Sawford, 1979: A model of katabatic winds. J. Atmos. Sci., 36, 619-630. Massman, W.J. and X . Lee, 2002: Eddy covariance flux corrections and uncertainties in long term studies of carbon and energy exchanges. Agric. For. Meteor., 113, 121-144. Meyers, T.P. and K.T. Paw U, 1986: Testing of a higher-order closure model for airflow within and above plant canopies. Boundary-Layer Meteorol, 37, 297-311. Miller, D.R., J.D. Lin and Z.N. Lu, 1991: Air flow across an alpine forest clearing: A model and field measurements, Agri. For. Meteor., 56, 209-225. Moeng, C.H., 1984: A large eddy simulation model for the study of planetary boundary-layer turbulence. J. Atmos. Sci., 41, 2052-2062. Morgenstern, K., T.A. Black, E.R. Humphreys, T.J. Griffis, G.B. Drewitt, T. Cai, Z. Nesic, D.L. Spittlehouse and N.J. Livingston, 2004: Sensitivity and uncertainty of the carbon balance of a  Pacific Northwest Douglas-fir forest during an El Nino/La Nina cycle. Agric. For. Meteorol. 123, 201-219. Novak, M.D., W. Chen, A.L. Orchansky and R. Ketler, 2000: Turbulent exchange processes within and above a straw mulch. I: mean wind speed and turbulent statistics. Agric. For. Meteorol., 102(2/3): 139-154. Ogura, Y. and N.A. Phillips, 1962: Scale analysis of deep and shallow convection in the atmosphere. J. Atmos. Sci., 19, 173-179. Oke, T.R., 1987: Boundary layer climate. 2 edition. Cambridge University Press, 435 pp. nd  Oliver, H.R., 1971: Wind profiles in and above a forest canopy. Q.J.R.M.S., 97, 548-553. Peltier, W.R. and T.L. Clark, 1979: The evolution and stability of finite amplitude mountain waves — II: Mountain wave drag and severe downslope windstorms. J. Atmos. Sci., 36, 1499-1529. Peltier, W.R. and T.L. Clark, 1983: Non-linear mountain waves in two and three spatial dimensions. Q.J.R.M.S., 109, 527-548. Plate, E.J., 1971: Aerodynamic characteristics of atmospheric boundary layers. AEC Crit. Rev. Ser., TID-25465. Raupach, M.R. and A.S. Thom, 1981: Turbulence in and above plant canopies. Ann. Rev. Fluid Mech., 13, 97-129. Raupach, M.R., P.A. Coppin and B.J. Legg, 1986: Experiments on scalar dispersion within a model plant canopy. Part I: the turbulence structure. Boundary-Layer Meteorol. 35, 21-52. Raupach, M.R., R.A. Antonia and S. Rajagopalan, 1991: Rough-wall turbulent boundary layers. Appl. Mech. Rev., 44, 1-25. Raupach, M.R. and J.J. Finnigan, 1997: The influence on meteorological variables and surfaceatmosphere interactions. J. Hydrology, 190, 182-213. Redelsperger J. and T.L. Clark, 1990: The initiation and horizontal scale selection of convection over gently sloping terrain. J. Atmos. Sci., 47, 516-541. Reich, P.B. and P. Bolstad, 2001: Productivity of evergreen and deciduous temperate forests. In: J. Roy, B. Saugier, H.A. Mooney: Terrestrial Global Productivity. Physiological Ecology. Academic Press, San Diego, 245-283. Ross, J., 1975: Radiative transfer in plant communities. Vegetation and the Atmosphere. Vol. 1, J.L. Monteith ed., Academic Press, New York, 13-55. Ruck, B. and E. Adams, 1991: Fluid mechanical aspects of the pollutant transport to coniferous trees. Boundary-layer Meteorol., 56, 163-195. Saigusa, N., T. Oikawa and S. Liu, 1998: Seasonal variations of the exchange of C 0 and H 0 between a grassland and the atmosphere: an experimental study. Agric. For. Meteorol, 89, 131-139. 2  2  Sakai, R.K., D.R. Fitzjarrald and K.E. Moore, 2001: Importance of low-frequency contributions to eddy fluxes observed over rough surfaces. J. Appl. Meteorol., 40, 2178-2192. Sasamori, T., 1972: A linear harmonic analysis of atmospheric motion with radiative dissipation. J. Meteorol. Soc. Japan, 50, 505-517. Sarmiento, J.L. andN. Gruber, 2002: Sinks for anthropogenic carbon. Phys. Today, 208, 30-36. Schmid, H.P., 1994: Source Areas for Scalars and Scalar Fluxes. Boundary-layer Meteorol., 67, 293-318. Schimel, D.S., 1995: Terrestrial ecosystems and the carbon cycle. Global Change Biology, 1,77-91. Schimel, D., T.G.F. Kittel, S. Running, R. Monson, A. Turnipseed and D. Anderson, 2002: Carbon sequestration studied in Western U.S. Mountains. EOS, 83, No. 40, 445-449. Schuman, F.G., 1962: Numerical experiments with the primitive equations. Meteor. Soc. Japan, 85107. Schumann, U., H . Thomas, H. Hartmut, S. Helmut and V. Hans, 1987: A mesoscale model for the simulation of turbulence, clouds and flow over mountains: Formulation and Validation examples. Beitr. Phys. Atmos. 60, 413-446. Seginer, I., P.J. Mulhearn, E.F. Bradley and J.J. Finnigan, 1976: Turbulence flow in a model plant canopy. Boundary-Layer Meteorol. 10, 423-453. Shaw, R.H., G.D. Hartog and H.H. Neumann, 1988: Influence on foliar density and thermal stability on profiles of Reynolds stress and turbulence intensity in a deciduous forest. Boundary-Layer Meteorol, 45, 391-409. Smagorinsky, J., 1963: General circulation experiments with the primitive equations. I. The basic experiment. Mon. Wea. Rev., 91, 99-164. Smolarkiewicz, P.K., 1983: A simple positive definite advection scheme with small implicit diffusion. Mon. Wea. Rev., I l l , 479-486. Smolarkiewicz, P.K., 1984: A fully multidimensional positive definite advection transport algorithm with small implicit diffusion. J. Comput. Phys., 54, 325-362. Song, B., J. Chen, P.V. Desanker, D.D. Reed, G.A. Bradshaw and J.F. Franklin, 1997: Modeling canopy structure and heterogeneity across scales: from crowns to canopy. Forest Ecol. Manage, 96, 217-229. Staebler, R . M . and D.R. Fitzjarrald, 2004: Observing subcanopy C 0 advection. Agri. For. Meteorol, 122, 139-156. 2  Stannard, D.I., J.H. Blanford, W.P. Kustas, W.D. Nichols, S.A. Amer, T.J. Schmugge and M.A. Weltz, 1994: Interprtation of surface flux measurements in heterogeneous terrain during the Monsoon's 90 experiment. Wat. Resor. Res., 30(5): 1227-1239.  Stephens, G., 1984: Short-wave parameterization revised to improve cloud absorption. J. Atmos. Sci., 41, 687-690. Stohl, A., S. Eckhardt, C. Forster, P. James, N . Spichtinger and P. Seibert, 2002: A replacement for simple back trajectory calculations in the interpretation of atmospheric trace substance measurements. Atmospheric Environment, 36, 4635-4648. Stull, R.B., 1988: An introduction to boundary layer meteorology. Kluwer Academic, 666 pp. Stull, R.B., 2000: Meteorology for Scientists and Engineers. 2nd Edition. Brooks/Cole, 502 pp. Sun, J., R. Desjardins, L. Mahrt and I. MacPherson, 1998: Transport of carbon dioxide, water vapor, and ozone by turbulence and local circulations. J. Geophys. Res., 103, 25873-2588. Taylor, P.A., P.J. Mason and E.F. Bradley, 1987: Boundary layer flow over low hills (a review). Boundary-Layer Meteorol.,39, 107-132. Tans, P., P. Bakwin and W. Guenther, 1996: A feasible Global Carbon Cycle Observing System: a plan to decipher today's carbon cycle based on observations. Global Change Biology, 2: 209-318. Thorn, A.S., 1971: Momentum absorption by vegetation. Q.J.R.M.S., 97, 414-428. Twine, T.E., W.P. Kustas, J.M. Norman, D.R. Cook, P.R. Houser, T.P. Meyers, J.H. Prueger, P.J. Staris and M.L. Wesely, 2000: Correcting eddy-covariance flux underestimates over a grassland. Agric. For. Meteorol., 103, 279-300. Uchijima, Z., 1961: On characteristics of heat balance of water layer under paddy plant. Bull. Nat. Inst. Agric. Sci., Vol. A, 243-263. Valentini, R., G. Matteucci, A.J. Dolman, E.D. Schulze, C. Rebmann, E.J. Moors, A . Granier, P. Gross, N.O. Jensen, K. Pilegaard, A. Lindroth, A. Grelle, C. Bernhofer, T. Grunwald, M . Aubinet, R. Ceulemans, A.S. Kowalski, T. Vesala, U . Rannik, P. Berbigier, D. Loustau, J. Guomundsson, H . Thorgeirsson, A.Tbrom, K. Morgenstern, R. Clement, J. Moncrieff, L. Montagnani, S. Minerbi and P.G. Jarvis, 2000: Respiration as the main determinant of carbon balance in Euopean forests. Nature, 404, 861-865. Wang, H . and E.S. Takle, 1996: On three-dimensionality of shelterbelt structure and its influences on shelter effects. Boundary-Layer Meteorol., 79, 83-108. Webb, E.K., G.I. Pearman and R. Leuning, 1980: Correction of flux measurements for density effects due to heat and water vapor transfer. Q.J.R.M.S., 106: 85-100. Wilson, N.R. and R.H. Shaw, 1977: A higher order closure model for canopy flow. J. Appl. Meteorol. 16, 1197-1205. Wilson, J.D., 1988: A second-order closure model for flow through vegetation. Boundary-Layer Meteorol. 42, 371-392. Wilson, J.D., J.J. Finnigan and M.R. Raupach, 1998: A first-order closure for disturbed plantcanopy flows, and its application to winds in a canopy on a ridge. Q.J.R.M.S. 124: 705-732.  Wofsy, S.C., M.L. Goulden, J.W., Munger, S.M. Fan, P.S. Bakwin, B.C. Daube, S.L. Bassow and F.A. Bazzaz, 1993: Net exchange of C 0 in a mid-latitude forest. Science, 260, 1314-1317. 2  Yamada, T, 1982: A numerical study of turbulent airflow in and above a forest canopy. J. Meteor. Soc. Japan, 60, 439-454. Yang, P.C., T.A. Black, H.H. Neumann, M.D. Novak and P.D. Blanken, 1999: Spatial and temporal variability of C 0 concentration and flux in a boreal aspen forest. J. Geo. Res. 104 (D22), 2765327661. 2  Y i , C , K.J. Davis, P.S. Bakwin, B.W. Berger and L.C. Marr, 2000: Influence of advection on measurements of the net ecosystem-atmosphere exchange of C 0 from a very tall tower. J. Geo.Res. 105, 9991-9999. 2  Zeng, P.T. and H. Takahasshi, 2000: A first-order closure model for the wind flow within and above vegetation canopies. Agric. For. Meteorol, 103, 301-313.  WFIS Model Dynamics A . l Momentum equations in spherical geometry The equations of motion in the geo-spherical system ^.-longitude, ^-latitude, r-radius that describe a balance between inertial, pressure- gradient, Coriolis, buoyancy and frictional effects are: Du Dt Dv  tgd l — w v + — uw = r r tgtj> 1 2  + — M  Dt Dw — Dt  +-VW =  r  1  r ,  2  2n  (K + V )= 2  2  r  1  dp , l (V-t), —+fv-f w + — r p r cos ip dX p \ \ dp . (V.r)2+-^r,,+-r23 —-fu + — r r p r dip p l  r  r> (J  +  B+f  g  1  y* u  +p  (V-r) --(r +r ) 3  (A-l) (A-2)  r  14'/' ^- + JL p dr p  ri2+-7l3  n  2 2  (A-3)  where velocity components are u (eastward), v (northward), and w (along the radius). B is buoyancy when phase changes exist, p is air density. (j = d(\np)/dz,f  = 2nsm0, a n d / * = 2Q. cos <p.  (A-4)  is the Reynolds stress tensor, which is parameterized using a first-order theory to represent sub-grid scale turbulence. The integration domain is bounded above by a rigid lid, which has the effect of reflecting vertically propagating waves. In order to suppress upperlid reflections, an artificial damping called Rayleigh friction is allowed to have finite values in the upper several levels.  A.2 Continuity equation The continuity equation adopts the 'anelastic' approximation (V • (p V) = 0), and is taken as: a  ^•{pu) + — — {cos <j) pv) + \—(r r cos (j> dX r cos^ dtp r dr 1  2  pw) = 0  (A-5)  A.3 Terrain-following coordinates Equation ( A - l ) through (A-5) represents the smooth earth equations. To convert to terrainfollowing coordinates, the vertical coordinate transformation is taken as Schumann et al. [1987]:  z = F{Q[\-hlH]  +h  (A-6)  where h = h(x,y) is the height of the topography, z is the Cartesian height, C, is the transformed vertical height in the new non-orthogonal coordinate system, and H is the height of the topmost fixed ' l i d ' of the model. F(Q is chosen as a monotonically increasing function of £ The boundary conditions on F(Q are taken as |  F ( 0 )  =°  (A-7)  In the case of F (C, ) = £, the above system reduces to the coordinates as used in Gal-Chen and Somerville [1975] and Clark [1977]. C, = 0 is the surface, £ = H defines the horizontal surface at the model top.  WFIS Model Physics B.l Conservation equation for potential temperature The conservation equation for 8 is taken as: — +— (u6) + — (v0) + — {w0) = S +— +—ot ox oy dz ox ay L  g  where S represents the source of heat, and H\, H g  + —~ oz  (B-l)  and Hj, represent sub-grid-scale  2  turbulence heat flux. ( B - l ) has been transformed to traditional meteorological coordinates from geospherical system in which it is convenient to formulate the second-order finite differences.  B.2 Radiation cooling schemes The code contains several options to calculate the infrared radiation under the absence of canopy. Here I only describe the R - C scheme [Redelsperger and Clark, 1990], which is the scheme used in this study. The general form of source-sink term in the 9 equation due to the long-wave radiative cooling processes is determined by the radiative flux divergence:  ( f ^ - ^ J ^ - ^ )  (B-2)  where F a n d F± are the total upward and downward radiation fluxes. C, is the height above r  the terrain.  B.2.1 Most detailed approach Using the broadband emittance form, the upward and downward long-wave radiative fluxes in (B-2) at height C, can be written respectively as [Stephens, 1984]: F (0 = or/ t  =0  + face)  y  daT  ]  dC  (B-3)  etc,  F & ) = - j ^ r , o ^ ° < e'V.O  = [A {Z,t> /^dv * doT {Q ) d B  v  (b-4) (B-5)  where A and B are the absorptivity and the Planck functions for wave number v, v  v  respectively. Substituting (B-3) and (B-4) into (B-2), one gets the model's most detailed procedure to calculate the long-wave radiation of broadband emittance: ( ^ ^ - ^ ^ { l o ^ O ^ ^ ^ W ^ O ^ ^ ^ ' }  (B-6)  In the troposphere, the principal effects of absorption and scattering are due to the rotational bands (18-100 urn) and vibration-rotational band (6.3 urn) of water vapor, the carbon dioxide band (around 15 urn), and the continuum absorption in the atmospheric window (8-14 um). These effects together with the cloud-water effect are included in the present model. Also,  e'(C,C)  is expressed as a function of path length for water vapor and  carbon dioxide, which are computed from integrals between C, and C,'. Therefore, the detailed approach takes the effect of carbon dioxide and water vapor as well as liquid water into account, gives a fairly accurate radiation-cooling budget. However, it requires computing double integrals that are obviously too expensive to use in a threedimensional dynamical model at every time step. A simpler and more economical formulation is the Sasamori radiation scheme.  B.2.2 Sasamori scheme The radiative effect at height C, may be decomposed into effects due to radiative exchanges between the height C, and (i) space, (ii) the ground, and (iii) other layers. Assuming that the latter effect is much smaller than the others, Sasamori (1972) proposed to neglect this term. Also assuming that the whole atmosphere has same temperature (isothermal atmosphere), except for infinitesimal layers near each boundary, leading to: it .&)R  =- ^—{o-(T\O-T m^- '(i,0)-aT (O^- U,^)} 4  7  •  (B-7)  4  £  (Tpcp)  £  %  .  <%  However, the Sasamori scheme (B-7) fails to give a sufficiently accurate estimation of the radiative flux divergence for the full troposphere due to the overly simplified approximation.  B.2.3 Combined approach (R-C scheme) The default approach used in this study is to combine the detailed approach with the Sasamori approach. This combined approach assumes that the radiative exchanges between different layers of the atmosphere cannot be neglected, but the effect of their temporal variations may be neglected to determine the radiative divergence. The combined radiative term is written as: r  a  d9\  SA  \a)  + R  \a J  (B-8)  R  On the right-hand side of (B-8), the first term is computed at each time step by (B-7). The term [ ] represents the difference between the complete form (detailed approach, (B-6)) and the approximate form (Sasamori scheme) (B-7). The detailed term is computed infrequently to keep the computational cost to a minimum. The combined scheme accounts for the effect of carbon dioxide and water vapor as well as liquid water, and also is computationally affordable. In the study, we use the combined scheme to calculate the infrared radiation in the heat budget of the model.  Sensitivity experiments and preliminary results C.l Sensitivity to tree-drag coefficient As mentioned in Chapter 2, it is difficult to measure the effective drag coefficient for a forest canopy, due to the shelter effect [Thom, 1971] and leaf orientation. The effective C„d of a single leaf measured in wind tunnel changes with the leaf orientation, turbulent scales, and intensity around the leaf [Raupach and Thom, 1981]. We cannot predict the drag on the foliage from the knowledge of the canopy geometry and the behavior of the plant elements in isolation. The following sensitivity experiments test the effect of a range of drag coefficients on the lee flow.  C . l . l Experimental design The tree-drag experiments employed a simple initialization scheme using an idealized dry atmospheric profile. Effects of radiation, moisture and surface heating were ignored. A background flow over a 1300 m high ridge, representing a simplified, smoothed cross section of the mountain to the west of Campbell River, was simulated. The atmosphere was taken as weakly stable, with Brunt-Vaisala frequency ./V = 0.001 s" from z = 0 to 2 km, N = 1  0.01 s" from z = 2 to 13.5 km and N = 0.02 s" above 13.5 km. A rather stiff breeze of 10 1  1  m/s from the west was prescribed at all heights. The Coriolis force was turned off. The experiments employed two nested domains. Table C . l outlines the general characteristics of these experiments. Ay and Az represent the grid size in the horizontal and vertical, respectively. NY and NZ are the number of grids in y and z directions, respectively. At denotes the time step, and C j is the tree-drag coefficient. Duration represents the model n  simulation time. A n isolated canopy of trees 33 m high was specified from y = 46 to 61 km in the lee of the mountain. The horizontal outer domain ranges from 0 to 87 km, and the nested domain was from 37.8 to 81 km. The leaf area density, a(z), was specified using (2-4).  Table C . l : Characteristics of sensitivity to drag-coefficient experiments. A l l experiments used two nests with the second nest starting from y = 37.8 km. Slashes separate values for the two nested domains. Experiments Ay(m) NY NZ At (s) Az j Cnd Duration (min) (m) D R A G 1A 150/75 578/578 140/42 6/3 4.0/1.333 0.001 60. D R A G IB 140/42 150/75 578/578 6/3 4.0/1.333 0.01 60. DRAG1C 140/42 2.0/0.6667 150/75 578/578 6/3 0.1 60. D R A G ID 150/75 578/578 6/3 140/42 1.0 1.5/0.50 60. m  DRAGOD  150/75  578/578  6/3  n  140/42  4.0/1.333  0.0  60.  C.1.2 Results The obstacle effects on the background flow in these experiments result from both the mountain and the isolated stand. Since the atmosphere is near neutral below 2 km (N = 0.001), the air near the surface will readily flow over the mountain and the wind might accelerate slightly to compensate for the thinned flow channel. The parameters chosen for these experiments is one where the dynamics is pretty much potential flow near the surface but with vertically propagating waves above the A B L that are mostly hydrostatic. Figure C . l shows the streamlines for the outer domain of experiments D R A G O D (C„</ = 0, no tree-drag force) and D R A G ID (C„d = 1, strong tree-drag force). A vertically propagating internal gravity wave is evident between 2 km and 13.5 km. The wave crests tilt upstream with increasing altitude, which is consistent with a vertical flux of energy. These waves have a vertical wavelength: A = 2nUIN = 6 km where N = 0.01. It is the 6 km vertical wavelength that can be clearly seen in the upper levels of Fig. C . l . Comparing the amplitude of the gravity waves in these two cases, D R A G 1 D has smaller amplitude waves, which indicates part of the momentum has been absorbed by the tree canopy in the surface layer. The influence of canopy drag on the flow is shown in the slight vertical deflection of the lowest streamline in Fig. C . l b . The obstacle effect of an isolated canopy is more evident in Fig. C.2. To clearly show the impacts of tree-drag on the surface flow, we sketch the topography above the graph frame and use z to represent height above local ground level ( A G L ) . For D R A G O D (C„d = 0.0) case, the canopy drag is zero and it does not affect the flow (Fig. C.2a), so the flow keeps its potential-flow properties. The little irregularities in the streamlines for D R A G O D are a result of some inaccuracies in the model M P I treatment. This has been corrected  subsequent to running these experiments. For D R A G 1 D (C„d =1.0), the canopy has a significant obstacle effect on the flow (Fig. C.2b). The strong canopy drag results in rapid near-surface transfer of horizontal momentum, deflecting the streamlines on the upwind side to well above the canopy that then return to the surface well to the lee of the canopy. Fig. C.2b shows a flow separation occurring immediately downstream of the isolated stand because the isolated canopy in the lee acts as an obstacle distorting the approaching flow. The primary cause of the flow separation is the loss of mean kinetic energy absorbed by the canopy, resulting in an adverse pressure gradient leading to flow separation.  Y  (KM I  Y  (KHI  Figure C . l : Streamlines for experiments DRAGOD (a, C„ = 0.0) and DRAG ID (b, C the outer domain. d  37.8  48.6  59.4  Y  (KM)  70.2  81.0  37.8  nd  48.6  59.4  70.2  = 1.0) for  81.0  Y (KHI  Figure C.2: Inner domain streamlines for experiments DRAGOD (a) and DRAG1D (b). The light rectangle area from 46 to 61 km represents the isolated canopy. Heights (z) are above local ground level.  Fig. C.3 shows the horizontal velocity fields for the nested domain for experiments D R A G 1 A , I B , 1C, and I D . As the drag coefficient increases from 0.001 to 1.0, the reverse motions occur further upstream for larger C„d (  =  0.1 and 1.0), indicating that flow  separation near the surface immediately downstream of the isolated canopy occurs earlier for larger tree drag. As simulation time progressing, a flow reversal just downstream of the canopy occurs for all cases. The reversal extends vertically well above the canopy, indicating a form-drag response; i.e., horizontal momentum exchange between the canopy and the free atmosphere. The larger the tree drag, the stronger the reverse motion. To the lee of the canopy the maximum reverse motion has a magnitude as large as 4.7 m/s for C„d = 1.0, and the affected height extends to 100 meters, around three times the canopy height. The physical explanation of the flow separation is the dramatic increase in drag resulting in almost complete kinetic-energy absorption. Within the canopy, the reverse flows extend only to a height of 20 m; i.e., they remain in the trunk space. The lack of tree drag in the trunk space could cause the reverse motion to penetrate. In the upper canopy, the strong tree drag suppresses any reverse motion. The larger the drag coefficient, the deeper the penetration of the reverse motions. This indicates the vertical coupling of momentum into the canopy could also depend on the canopy density, which determines the canopy drag. Figure C.4 shows the vertical velocity for tree-drag experiments in the inner domain. The case of C„d = 0.001 is the closest to the pure mountain flow, with updrafts on the upwind side of the mountain and downdrafts in the lee. As C„d is increased, the contribution of the adverse flow resulting from the canopy drag shifts the region of downdrafts to upstream. The larger the drag, the further the shifting. Right at the downstream edge of the isolated canopy, weak upward motions indicate the presence of a weak reversed rotor for larger tree drag cases. We also notice little difference between C„d = 0.1 and C„d = 1.0. Because the time scale of flow adjustment within the canopy is already very small for C„d =0.1, further increasing the drag coefficient makes little difference for the resolution of the model. This suggests that a value of C„d between 0.1 and 1.0 is reasonable for our further numerical simulations. In summary, the tests of sensitivity to drag coefficients indicate that an isolated stand acts as an obstacle that deflects the approaching lee flow, causing upstream flow separation  and downstream flow reversal. The bigger the tree-drag, the earlier the flow separation and the stronger the reverse motion. A possible application of these results is that the reverse eddy motions might be expected for flow regimes such as behind windbreaks, near clearings or forest edges, or where the vegetation is sparse and irregular. Further, although there are some arguments on the effective drag coefficient, the tree-drag experiments suggest that the drag coefficient between 0.1 and 1.0 makes only minor difference to the disturbed flows.  37.8  4 8 . 6  Y  5 9 . 4 (KM J  7 0 . 2  8 1 . 0  3  7  - 8  4  8 - &  Y  5 9 . 4 (KM!  7 0 . 2  8 1 . 0  Figure C.3: Horizontal velocity for inner nest for experiment DRAG1A (C„d = 0.001), lB(C„d = 0.01), 1C (C„d = 0.1), and ID (C„</ = 1.0), respectively. Contour interval is 1.5 m/s. Solid contours represent flow from west to east. Dashed contours denote the reverse motions.  9 9  81  .0  81  99  .1  79  .3  59  .5  39  .6  C„ = d  1 9 . B  37  .0  1.0  1 9 .  .8  37 .  4 8 . 6  5 9 . 4  Y  7 0 . 2  81  (KM)  Figure C.4: Vertical velocity for the inner domain for DRAG1 [A, B, C, D]. Solid contours represent upward, dashed contours denote downward. Contour interval is 0.25 m/s.  C.2 Sensitivity to topography and upstream flow In this section, we perform another two series of experiments with various heights of a hill and with changed upstream flows. These tests represent a small sample of possible largescale flows. The purposes of these tests are to investigate the conditions under which we might expect tree-drag to induce downstream backward motions, i.e. flow reversal. Since horizontal variability of forests is natural and is becoming increasingly common due to forest-management  practices, we can expect forest wind breaks to play a role in local  dynamic flow behavior. Further, the variability of canopy drag and terrain as well as the resulting dynamics are also likely to play an important role on the interpretation of the eddy- covariance measurements and N E E estimates.  C.2.1 Experimental design The preceding tree-drag experiments showed there were only minor differences to the disturbed flow for C„d ranging between 0.1 and 1.0. With this in mind, we choose C„d = 0.15 for the following eight experiments outlined in Table C.2. This value was also suggested by Amiro [1990], and is considered typical for mixed deciduous forests. The topography experiments (HMNT) have the same shape of topography as in the tree-drag experiments except the maximum height of the mountain is chosen as 100 m, 500 m, 900 m and 1300 m. The upstream wind speed VQ is set to 10 m/s for the H M N T tests. The upstream experiments ( R G M V ) fix the height of the mountain to 1300 m but vary the upstream wind speed VQ from 1 m/s, 3 m/s, 5 m/s to 10 m/s, respectively. The atmospheric stratification condition is taken as the same as for the tree-drag experiments. The sensitivity experiments employ triple nested domains with horizontal grid sizes of 150 m, 75 m and 25 m, respectively. The minimum vertical grid spacing is 12 m, 6 m, and 3 m, respectively. The second nest starts fromy = 37.8 km to 81 km, the third domain starts from 43.2 km to 64.8 km. The same initial conditions as for the tree drag tests are used. Table C.2 outlines the characteristics of topography and upstream experiments. Ay and Az represent the grid size of horizontal and vertical, respectively. A^ represents the time step, VQ is the upstream wind speed. Table C.2: Characteristics of sensitivity to topography (HMNT) and to upstream flow (RGMV) experiments. Slashes separate values for the three nested grids. • ,4y(m) Az(m) Tests At (s) K (m/s) Height of hill (m) 0  HMNT1 HMNT2 HMNT3  150/75/25 150/75/25 150/75/25  12/6/3 12/6/3 12/6/3  7.5/3.75/1.25 6.0/3.0/1.0 4.0/2.0/0.67  10 10 10  100 500 900  HMNT4 (RGMV4) RGMV1  150/75/25  12/6/3  4.0/2.0/0.67  10  1300  150/75/25  12/6/3  8.6/4.3/1.4  1  1300  RGMV2  150/75/25  12/6/3  6.0/3.0/1.0  3  1300  RGMV3  150/75/25  12/6/3  7.5/3.75/1.25  5  1300  C.2.2 Results The features of canopy flow in the lee of mountains will be a result of both mountain forcing and local tree drag. These experiments, outlined in Table C.2, are still dominated by essentially potential flow near the surface that exhibits mountain-wave generation above the  A B L . We chose a range of mountain heights from nearly flat to steeply sloped, to see i f we could observe any asymptotic behaviors of the canopy flow. The effects of different mountain heights on the canopy flow are shown in Fig. C.5. The mean wind speed is significantly increased at the upwind side of the stand. When the air flows over the mountain, the increasing mountain height forces air through a narrower channel causing increased wind. The higher the mountain, the stronger the acceleration. Stronger winds approaching the isolated canopy result in stronger tree drag, increasing the flow deflection and causing earlier flow separation. A reverse eddy forms just downstream of the canopy. We also see from the topography experiments that reversed flow sets in when v  max  reaches about 16 m/s for H M N T 3 (Fig. C.5). The smoothness of the contours in H M N T 3 as opposed to H M N T 2 required we increase the eddy mixing to a value about 10 m /s for 2  these experiments. Otherwise the strong shears causes high-frequency oscillations, resulting from poor resolution of the local shear gradients. Upstream flow experiments R G M V [1-4] vary the upstream wind speed from 1,3, and 5 to 10 m/s, respectively. Figure C.6 shows the horizontal velocity for the second nest. To highlight the effects of mountain-wave flow, the contour interval in the plots is chosen to be proportional to the upstream wind speed, and is 0.1, 0.3, 0.5 and 1.0 m/s, respectively. We see a systematic increase in 'normalized' reverse flow as Vo is increased. If the entire flow were neutral we would expect to see very little change as Vo is increased. This increase suggests an increased normalized contribution from gravity-wave effects. In summary, the increasing mountain height forces air through a narrower channel, accelerating the approaching wind at the upwind side of the isolated stand. Stronger approaching winds results in stronger tree-drag and earlier flow separation and reversals. Increasing the upstream winds is also more favorable for the reversed motions. The interaction between the mountain and the canopy makes the properties of the flow more interesting and also more difficult to study. To interpret measurements made near hills and windbreaks one must consider all these factors, and distinguish which factor is dominant before reasonable explanation can be made.  43 . 2  48.6  54 . 0 T < KM I  59 . 4  Figure C.5: Horizontal velocity for the third nested domain for experiments HMNT1 (a, 100 m) HMNT2 (b, 500 m), HMNT3 (c, 900 m) and HMNT4 (d, 1300 m). Solid contours represent the direction of basic flow; dashed contours mark the reverse motion. Contour interval is 1.5 m/s.  Figure C.6: Horizontal velocity for the second domain for R G M V [1-4] experiments. Solid contours represent the direction of upstream flow. Dashed contours mark reverse motion. Contour interval is 0.1, 0.3, 0.5 and 1.0 m/s, respectively.  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

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

Comment

Related Items