Note

Please report issues or ask questions about this site on the GitHub page.

Boundary Layer Meteorology

Note

Please report issues or ask questions about this site on the GitHub page.

Glossary

AmeriFlux

The AmeriFlux network is a community of sites and scientists measuring ecosystem carbon, water, and energy fluxes across the Americas, and committed to producing and sharing high quality eddy covariance data.

EC

Eddy Covariance

FAO

Food and Agriculture Organisation, a specialised agency of the United Nations that leads international efforts to defeat hunger.

Fetch

Area influencing a sensor. Typically the area upwind of the sensor.

LAI

Leaf Area Index

Obukhov Length (\(L\))

A parameter with dimension of length that gives a relation between parameters characterizing dynamic, thermal, and buoyant processes. More detailed explanation refers to AMS wiki.

OHM

Objective Hysteresis Model

SEB

Surface Energy Balance

URAO

University of Reading Atmospheric Observatory

Note

Please report issues or ask questions about this site on the GitHub page.

Surface energy balance

The surface energy balance (SEB) can be defined for an extensive simple surface as (units: \(\mathrm{W \ m^{-2}}\)):

()\[Q^*= Q_H + Q_E +Q_G\]

where \(Q^*\) is the net all wave radiation, the turbulent sensible heat flux (\(Q_H\)) and the turbulent latent heat flux (\(Q_E\)) and the soil heat flux (\(Q_G\)) (units: \(\mathrm{W \ m^{-2}}\)).

Note other notation that is used for the SEB are (in the same order as above):

\[R_n= H + LE + G\]

Bowen ratio

The Bowen ratio (\(\beta\)) is

()\[\beta= Q_H / Q_E\]

Radiation balance

The \(Q^*\) (or \(R_n\)) within the SEB Surface energy balance consists of:

()\[Q^*= K_\downarrow- K_{\uparrow} + L_\downarrow- L_{\uparrow}\]

which includes the incoming (\(\downarrow\)) and outgoing (\(\uparrow\)) shortwave (K) and longwave radiation (L) fluxes.

SEB and Radiation Balance Measurements

Each of these fluxes can be directly measured or modelled using several methods (and data inputs).

Examples of instruments in the URAO are listed. From papers you will be able to determine how the fluxes and other variables are measured at the site you are studying.

Radiation fluxes

Various types of radiometers are used. For shortwave radiation, pyranometers are used, and for longwave radiation, pyrgeometers are used. The source area or field of view of a radiation sensor is fixed by geomtery.

Soil Heat Fluxes

Soil heat flux plates are used with temperature sensors above to determine the heat gain/loss between the plate (e.g. at 0.05 m below the surface) and the surface. In more complex environments the storage heat flux (heating and cooling) of the whole volume needs to be accounted for. For example in a forest, the trees (trunk, branches, leaves, air) as well as the soil itself. So in most environment the soil heat flux is one part of the storage heat flux.

Turbulent heat Fluxes

The turbulent heat fluxes and momentum can be measured using Eddy Covariance techniques (see Eddy Covariance). These require a sonic anemometer and an infrared gas analyser.

As the source area of EC sensors varies with wind direction, fetch, stability etc, the surface characteristics may change with time if the fetch is not homogeneous.

SEB and Radiation Balance Modelling

Generally, to calculate a convective or conductive flux, data are needed to determine the size of the gradient (e.g. temperature difference) and the ability of the medium (e.g. air, soil) to transfer heat (or mass). The latter may use a resistance scheme, an eddy diffusivity, or other approach, which changes with the state of the medium (e.g. stability if air, moisture state if soil).

One method to model the latent heat flux uses the Penman Monteith equation (Penman).

Note

Please report issues or ask questions about this site on the GitHub page.

Penman

The combination equation (Penman, 1948)

Given the difficulty in directly observing evaporative fluxes – and the need to estimate evaporation to determine irrigation of crops Penman (1948) developed the combination equation. Only standard meteorological observations at some height above either water or a surface covered with well-watered short grass, were required.

This equation “combines” two physical controls on evaporation: the supply of available energy \(Q^*-Q_G\) and the turbulent diffusion of water vapour from the surface. Firstly, the sensible and latent heat fluxes should be written in resistance notation

()\[ Q_{H} = - \rho c_{p}\frac{T - T_{0}}{r_{h}}\]

()\[ L_{v}E = - \frac{L_{v}}{R_{v}T}\frac{e - e_{0}}{r_{v}}\]

We can rewrite Eq.6.2 in terms of \(r_h\) and the effective psychrometer constant \(\gamma^*\):

()\[L_{v}E = - \rho c_{p}\frac{e - e_{0}}{{\gamma^*}r_{h}},\]

where the “thermodynamic” psychrometer constant is given by \(\gamma = \rho c_{p}/L_{v}\varepsilon\) and has a value of 0.66 hPa K^-1 near sea level, the ratio of gas constants is \(\varepsilon=R_a/R_v\), and the effective psychrometric constant is \(\gamma^*=\gamma r_v/r_h\).

The key to estimating the exchange of sensible and latent heat with the air above a surface is to establish the surface temperature \(T_0\). If this is not known, then it can be eliminated from the equations by assuming that the saturation vapour pressure is a linear function of temperature for small temperature differences. Additionally, it is assumed that the surface vapour pressure is at saturation, i.e.

()\[e_{0}=e_{s}\left(T_{0}\right)=e_{s}(T)-s\left(T-T_{0}\right)\]

By substitution of Eq.6.4 into Eq.6.2 we get \(Q_E\) in terms of \(T-T_0\) which is not known. But by combining the new equation for \(Q_E\), \(T-T_0\) may be eliminated from Eq.6.1 for \(Q_H\).

The final assumption is that of energy balance closure (EBC), i.e. that \(Q^*-Q_G=Q_H+Q_E\). By substitution of the new equation for \(Q_H\), we obtain the Penman equation

()\[L_{v}E = \frac {s\left( Q^* - Q_{G} \right) + \rho c_{p}\left( e_{s}(T) - e \right)/r_{H}} {s + \gamma_{*}}\]

giving latent heat flux as a function of variables measured at one height only.

Resistance notation

In the combination equation, fluxes are recast in terms of aerodynamic resistances. Thus, the momentum flux across a finite layer \(\Delta z = \left( z_{2} - z_{1} \right)\) can be expressed as

()\[\tau = \frac{\rho\left( u_{2} - u_{1} \right)}{r_{a}}\]

where \(r_a\) is the aerodynamic resistance. If \(u_1\) is taken to be the \(u(z_0)=0\), then the resistance can be related to the log law

()\[r_{a} = \frac{u(z)}{u_{*}^{2}} = \frac{ln\lbrack\frac{z - d}{z_{0}}\rbrack}{ku_{*}} = \frac{\left\lbrack ln\lbrack\frac{z - d}{z_{0}}\rbrack \right\rbrack^{2}}{k^{2}u(z)}\]

Notice that the quantity in square brackets depends only on the nature of the surface: if that is known, then \(r_{a}\) is determined only by the single measurement of mean wind speed. For fluxes this leads to an expression predicting a flux of some variable \(\alpha\), \(F_{\alpha}\) in the form

()\[F_{\alpha} = u_{*}^{2} \left( \frac{\phi_{m}}{\phi_{\alpha}} \right) \frac{\Delta\alpha}{\Delta u} \equiv \frac{\Delta\alpha}{r_{\alpha}},\]

in which \(r_{\alpha} = \left( \frac{\Delta u}{u_{*}^{2}} \right) \times \left( \frac{\phi_{\alpha}}{\phi_{m}} \right)\) is the aerodynamic resistance of the layer (some authors use conductance, which is the reciprocal of \(r_{\alpha}\)). In the case of momentum, \(\phi_{\alpha}=\phi_m\). For non-neutral conditions:

()\[r_{a} = \frac{\left\lbrack \ln\left( \frac{z - d}{z_{0}} \right) - \ \psi_{m}\ \left( z^{'}/L \right) \right\rbrack^{2}} {k^{2}u}.\]

Under stable conditions, the aerodynamic resistance for sensible heat transfer (\(r_h\)) is usually taken as equal to \(r_a\). Under unstable conditions, the assumption \(\phi_{h}=\phi_{m}^{2}\) means that we need to multiply \(r_a\) by \(\phi_{m}\) to get \(r_h\). In both cases, the aerodynamic resistance for water vapour and other scalar fluxes is generally assumed to be the same as that for sensible heat. Note that, in general, \(z\) should be replaced by \(z-d\) if the sampling height is not much greater than the zero-plane displacement height \(d\).

The Penman-Monteith equation

Evaporation from a vegetated surface is also impacted by plant physiology. Stomata in the leaves open to allow transfer of CO₂ during photosynthesis, but can close when the plant is water stressed to avoid undue loss of water vapour. Monteith (1965) adapted Penman’s equation to allow for this effect, giving what we now know as the Penman-Monteith equation where the symbols have the same meaning, except that \(\gamma^*\) is an apparent psychrometer ‘constant’. Over crops the resistance to evaporation is larger than the resistance to heat transfer, due to canopy resistance. To allow for this, the effective psychrometer constant is usually assumed to be of the form

()\[\gamma_{*} = \gamma\left( \frac{r_{h} + r_{s}}{r_{h}} \right) = \gamma\left( 1 + \frac{r_{s}}{r_{h}} \right),\]

where \(r_s\) is an effective surface resistance. The latter depends in a complicated way on soil moisture, type of vegetation, fractional cover, and time of year. It is usual to consider it as the result of a canopy (or crop) resistance (\(r_{sc}\)) and a bare soil resistance (\(r_{ss}\)) acting in parallel, so that

()\[\frac{1}{r_{s}} = \frac{\left( 1 - A \right)}{r_{\text{sc}}} + \frac{A}{r_{\text{ss}}}\]

where \(A\) is an effective fraction of bare soil area. Appropriate values of crop resistance are known for various types of vegetation. For ‘moist’ surface conditions during the day, \(r_{ss}\) is usually taken to be 100 s m^-1.

Finally, \(r_h\) can be approximated by the aerodynamic resistance \(r_a\) as it is more readily measurable. For observations close to the ground (e.g. below 3 m) the stability correction can be neglected, and \(r_a\) therefore estimated using (4). The Penman approach does not require ‘special’ equipment but assumptions about the transfer processes at the Earth’s surface as well as the nature of the surface itself are needed.

The P-M equation is used practically by the FAO. Evapotranspiration is computed for a grass reference crop (\(ET_0\)), then multiplied by a crop coefficient (\(K_c\)) to give an estimate of actual evapotranspiration. For instance, \(K_c \sim 1 \textrm{ to } 1.2\) for cabbage, but 1.1 to 1.5 for sugar cane. To aid calculation of irrigation requirements, \(ET_0\) is expressed in \(\textrm{mm day}^{-1}\) and can range from \(1–3\textrm{ mm day}^{-1}\) in cool, arid regions to \(5–6\textrm{ mm day}^{-1}\) in warm tropical regions.

Penman Monteith Method Measurements

Meteorological variables measured at one height can be used to estimate the evaporative flux from a surface using the Monteith (1965) adaptation for vegetated surface of Penman (1948):

()\[Q_{E} = \frac{s\left( Q^{*} - Q_{G} \right) + \rho c_{p}(e_{s} - e)/r_{H}} {s + \gamma\ (1 + r_{s}/r_{a})}\]

For the practical, some assumptions are made in addition to those made in deriving the equation. The aerodynamic resistance for heat transfer is assumed to be equal to that for momentum transfer, i.e. \(r_h \sim r_a\).

Aerodynamic Resistance

To calculate \(r_a\), one could assume neutral conditions are applicable or apply stability correction. Surface characteristics influence the surface (\(r_s\)) or canopy resistance \(r_c\). Around the site the surface characteristics vary with grass in the near fetch but arrange of other land covers further away from the sensors.

Canopy or surface resistance

By re-arranging the PM equation, with EC and other observations you can determine the surface resistance (\(r_s\)) or its inverse surface conductance (\(g_s\)) (Ward et al. 2016):

()\[\frac{1}{g_{s}} = r_{s} = \left( \frac{s}{\gamma}\frac{Q_{H}}{Q_{E}} - 1\ \right) r_{\text{av}} + \frac{\rho c_{p}\text{VPD}}{\gamma Q_{E}}\]

where \(\text{VPD}\) (\(=e_s-e\)) is the vapour pressure deficit. For our purposes we will assume \(r_h \sim r_a \sim r_{av}\).

Note

Please report issues or ask questions about this site on the GitHub page.

Model Parameters

Land surface models use parameters to describe the surface. For example to model the latent heat flux using the Penman Monteith equation the following parameters are needed.

Albedo

From the short-wave radiation (K), within the Eq.3.1 Radiation balance the albedo (\(\alpha\)) is calculated:

()\[\alpha= K_{\uparrow} / K_\downarrow\]

using the incoming (\(\downarrow\)) and outgoing (\(\uparrow\)) shortwave radiation (K) fluxes.

Typical Values

Note

The table below is based crowd-sourced dataset at Albedo Collection. Please report issues there if any found.

Calculated values of albedo at selected AMF sites

site	location	land cover type	day of year	time of day	value	reference	remarks
US-SRG	Santa Rita Grassland	Grassland	27	12:00	0.24	Moore (1976)	summer value of grassland
US-VAr	Vaira Ranch-Iona	grasslands	10	12:00	0.45	Ryu et al. (2008)
US-MMS	Morgan Monroe State Forest	Deciduous Broadleaf Forest			0.13	Liu (2017)	value has seasonal variability-approximated from summer values
Crop field	UK	Cropland	25	13:00	0.2	Page (2012)	a short growing crop
CA-Qsu	Quebec-Eastern Boreal	Black Spruce / Jack Pine Cutover	15th December	12:00	0.907	Betts and Ball (1997)	winter time
CA-Qfo	Quebec	Eastern Boreal - Mature Black Spruce	196	12:00	0.081	Betts and Ball (1997)	none
40-50N	40-50N	Deciduous Broadleaf Forest	74	12:00	0.142	Gao et al. (2005)
Saskatchewan	Western Boreal	Mature Black Spruce	annual mean	daily mean	0.145	Rosbjerg (1997)
US-MOz	Missouri USA	Deciduous Broadleaf Forest		Midday	0.15+/-0.02	Moore et al. (1996)	Valid for summer days (a case study of a Deciduous forest)
US-Oho	Oak Openings USA	deciduous broadleaf forest	July		0.159	Strugnell et al. (2001)
CA-Obs	Saskatchewan	ENF			0.260	Bright et al. (2014)	January average

Roughness length (\(z_0\)) and displacement height (\(d\))

If the displacement height is known, or is negligible, the logarithmic law equation can be rearranged with observed \(z_0\) and mean wind speed to allow \(z_0\) to be determined. As this may vary we normally take median of a minimum of 20 results for a wind direction sector. If you have a period with a lot of neutral conditions you may be able to get a lot of samples rapidly.

\[𝑧_0 = (𝑧−𝑑) exp ⁡[−(𝑈_𝑧 𝜅)/𝑢_∗ ]\]

Literature values of roughness parameters collected by students of BLM class

group	land cover	z0	zd	ra	reference
1-1
1-2
2-1
2-2
3-1	Grasslands		0.25m	10-20	de Miguel and Bilbao (1999)
3-2	Suburban Neighbourhood Park	0.03	0.2		Kent et al. (2017)
4-1
4-2	Boreal, Mature Black Spruce	0.22	9.66	50	Kettridge et al. (2013)
5-1
5-2	Broadleaf Deciduous Forest	2.9	20.1		Maurer et al. (2015)
6-1
6-2	Cropland	0.062		77.0	Liu et al. (2007)
7-1
7-2
8-1
8-2	Grasslands	0.026	0.1	40	Dunin et al. (1978)
9-1
9-2
10-1
10-2

Calculated values of roughness parameters at selected AMF sites

site	z0	zd	ra	site description
CA-Obs
US-Blk
US-MMS
US-MOz	1.48	15.40	9.02	Deciduous Broadleaf Forest
US-Dia	0.04	0.65	19.1	Grassland
US-KUT	0.03	0.049	159.6	Temperature grassland
CA-Qcu	0.24	8.4	27.7	Boreal Black Spruce/Jack Pine Cutover
CA-Qfo	1.86	9.66	8.55	Boreal Mature Black Spruce
US-Slt
US-UMB	2.32	15.4	3.4	Deciduous Broadlead Forest
US-Br3
US-Bo1	0.04	1.37	50.1	Croplands
US-NC1	0.38	3.5	23.5	Young pine forest
US-Whs
US-SRG
US-Var	0.09	0.03499	58.0899	Grasslands
US-Bo2
US-Br1
CA-TPD
US-Oho	2.66	1.75	5.63	Deciduous broadleaf forest

How does it vary with wind direction?

A rule of thumb for calculating d is to assume it is ~0.7 \(h\) where \(h\) is the height of the canopy. As the heights may vary with direction you can determine how much this may vary. What are expected to be consistent sectors?

The wind profile can also be used to determine \(z_0\) and \(d\) if there are more than 2 levels in the profile. This requires fitting a straight line (linear regression) through the data to determine the intercept, which provides the \(z_0+d\) value. See equations 1-2 in [9].

For References see list

Note

Please report issues or ask questions about this site on the GitHub page.

LAI (Leaf Area Index)

Key Facts

• Varies with phenology (state of the vegetation e.g. leaf on or leaf off, bud burst).

• Measure of the area of leaf surface relative to the ground area. If you look up at an adult, healthy tree there will be typically many layers of leaves between the ground and the top of the canopy. Hence the LAI values are typically greater than 1 m² m^-2.

• Units: m² m^-2

Porosity

If air can blow through an object (e.g. a tree) it is more porous than those one the that air can not (e.g. buildings). The latter are referred to as bluff bodies. Trees porosity may change through the year with phenology. For example, when deciduous plants lose their leaves (fall off) they become more porous.

Measuring LAI

Direct measurements

• Destructive sampling - measuring the area of all the leaves of a tree

• Collect all the leaves that fall off and measuring their area

Indirect measurements

• Measurements (above and below the canopy) of radiation (e.g. using LAI-2000, LI-COR) which are then used to infer LAI

• Satellite based products e.g. MODIS)

Model LAI

• Depends on:

  • temperature (growing degree days)

  • solar radiation/day length

Parameters LAI influences

• Albedo

• Roughness length and displacement

• Aerodynamic resistance

• Surface/ canopy resistance

For assignment we needs a simple measure of the value and how it varies. For example, do you expect the LAI to vary through the year for the land cover you are looking at? If yes, when should it increase/ decrease? What is an approximate maximum (typical value)?

Note

Please report issues or ask questions about this site on the GitHub page.

Stability

Modifications to the logarithmic profile are required in conditions of non-neutral stability, using the results of Monin-Obukhov theory. This theory of the surface layer derives relations between the vertical variation of wind speed u(z) and potential temperature \(\theta(z)\) (which approximates the measured temperature T close to the surface), the scaling factors for momentum and temperature, \(u^*\) and \(T^*\), and the Monin‑Obukhov stability parameter

()\[\frac{z'}{L} = - \frac{k\left( z - d \right)\frac{g}{\theta_{0}}\frac{H}{\rho c_{p}}}{u_{*}^{3}}\]

where \(L\) is the Obukhov length and \(z^{’}= z - d\). NB: the surface temperature \(\theta_0\) is an absolute temperature (units: K). The logarithmic profile relation can be rewritten for wind speed to include the stability corrections

()\[\bar{u}(z) = \frac{u_{*}}{k}\left\lbrack \ln\left( \frac{z - d}{z_{0}} \right) - \Psi_{m}\left( \frac{z - d}{L} \right) + \Psi_{m}\left( \frac{z_{0}}{L} \right) \right\rbrack\]

and similarly, for potential temperature:

()\[\bar{\theta}(z) = \theta_{0} + \frac{T_{*}}{k}\left\lbrack \ln\left( \frac{z - d}{z_{h}} \right) - \Psi_{h}\left( \frac{z - d}{L} \right) + \Psi_{h}\left( \frac{z_{h}}{L} \right) \right\rbrack\]

where the turbulent temperature scale \(T_*\) is given by \(T_{*} = - \overline{w^{'}T^{'}}/u_{*} = - Q_{H}/(\rho C_{p}u_{*})\), \(\Psi_{m}\) is the integral stability correction function for momentum and \(\Psi_{h}\) is the integral stability correction function for heat.

There are a number of forms of \(\Psi_{m}\) and \(\Psi_{h}\); one set of forms from Foken (2008) are as follows:

under unstable conditions:

()\[\begin{split}\begin{array}{c} {\psi_{m}(\zeta)=\ln \left[\left(\frac{1+x^{2}}{2}\right)\left(\frac{1+x}{2}\right)^{2}\right]-2 \tan ^{-1} x+\frac{\pi}{2} \text { for } \frac{z}{L}<0} \\ {\psi_{h}(\zeta)=2 \ln \left(\frac{1+y}{2}\right) \text { for } \frac{z}{L}<0} \end{array}\end{split}\]

with \(x=(1-19.3 \zeta)^{1 / 4} \quad y=0.95(1-11.6 \zeta)^{1 / 2}\).

under stable conditions:

()\[\begin{split}\begin{array}{l} {\psi_{m}(\zeta)=-6 \frac{z}{L} \quad \text { for } \quad \frac{z}{L} \geq 0} \\ {\psi_{h}(\zeta)=-7.8 \frac{z}{L} \quad \text { for } \quad \frac{z}{L} \geq 0}\end{array}\end{split}\]

Note that both \(T_*\) and \(z’ / L\) have the opposite sign to \(Q_H\) (which is positive in unstable conditions and negative in stable conditions). If \(z’/z_0 \gg 1\) then the third term can assumed to be negligible (Garratt 1992).

Other stability metrics include the Richardson number:

• Gradient

• Bulk

• Flux

Bulk Richardson number is the ratio of thermally produced turbulence and turbulence generated by vertical shear or the ratio of free or forced convection (thermal: mechanical)

()\[R_{B}=\frac{\left(g / T_{v}\right) \Delta \theta_{v} \Delta z}{(\Delta U)^{2}+(\Delta V)^{2}}\]

where \(g\) acceleration due to gravity, \(T_V\) virtual temperature, \(\Delta \theta_{v}\) change (difference) in potential temperature, \(\Delta z\) change in height \(\Delta U\) change in \(U\) wind-speed, and \(\Delta V\) change in \(V\) wind-speed.

Note

Please report issues or ask questions about this site on the GitHub page.

Eddy Covariance

Sonic Anemometer

A sonic anemometer derives the wind-speed from transit times of acoustic pulses travelling in both directions along a fixed path. The wind-speed component along the path is proportional to the difference between the transit times. Three sets of transducer pairs are used to derive the three components of the wind vector u, v, w. In addition, the speed of sound can be deduced:

()\[c^2=403T_s,\]

where \(T_s\) is the sonic temperature:

()\[T_s=T(1+0.32 e/p)\]

where \(T\) is the dry bulb temperature, \(e\) is the water vapour pressure, and \(p\) is the air pressure. The sonic temperature is approximately equal to the virtual temperature \(T_v\).

In practical use, the covariance: \(\overline{w^{\prime} T_{s}^{\prime}} \approx \overline{w^{\prime} T^{\prime}}\) and therefore the sonic anemometer can be used to estimate turbulent sensible heat flux (\(Q_H\)). Sonic anemometers are the standard instrument used to observe turbulence in the atmospheric boundary layer.

Principle behind sonic measurement of velocity and speed of sound

Infra-red Gas Analyser (IRGA)

IRGA allows measurement of variations in water vapour and carbon dioxide allowing the latent heat flux and carbon dioxide to be measured. The specific humidity of water vapour \(q\) is expressed in units of kg kg^-1. The absolute humidity (kg m^-3) is derived by taking the molecular weight of water into account (1 mol = 18 g = 0.018 kg) and similarly for carbon dioxide concentrations (molar mass 44 g mol^-1).

1. Open-path IRGA: Measurements are made at station pressure.

2. Closed-path IRGA: Air is sucked down a tube into the instrument itself

Schematic of the LI-COR Li 7500 open path infra-red gas analyser (IRGA) (source: LI-COR)

Instruments at Different Sites

Gill R3 sonic anemometer (source: Gill Instruments).

1. Site: URAO

   • see Observatory.

2. Site: AmeriFlux

   • Each site has a key reference that gives the details of the instruments used.

   • Please see AmeriFlux Site page for details (login required with free registration).

Eddy covariance (EC) method to measure fluxes

Note

\(\tau\) unfortunately (by convention) is used for both momentum and time lag

EC is regarded as the standard method to measure fluxes as it quantifies directly the contribution by each eddy. EC fluxes of momentum, heat and moisture require fast response (10-20 Hz) measurements to capture the scales of turbulence contributing to fluxes. The momentum flux (or Reynolds stress) is calculated from:

()\[\tau=\rho \sqrt{\left(\overline{u^{\prime} w^{\prime}}^2 +\overline{v^{\prime} w^{\prime}}^2\right)}\]

This considers any vertical transfer of lateral momentum \(\overline{u^{\prime} w^{\prime}}\). The sensible heat flux is given by:

()\[Q_{H}=\rho c_{p} \overline{w^{\prime} T^{\prime}}\]

assuming that \(T\approx \theta\) , and the latent heat flux (\(Q_E\) ) is given by:

()\[Q_{E}=L_{V} \overline{w^{\prime}} q^{\prime}\]

where \(L_{V}\) is the latent heat of vaporisation (\(\approx 2.45 \times 10^6 \text{kg}^{-1}\)) and q is the specific humidity. These equations assume the vertical component of the flux dominates, i.e. flow is homogeneous and steady.

Rotation of EC data

In theory, the wind component u is defined as being aligned with the mean wind direction, and thus the mean vertical component. The instrument itself has a fixed frame of reference, so how is this achieved? The frame of reference is rotated to align the new u axis with the measured mean wind vector.

This is called double rotation as it is usually done in two steps:

1. rotate through angle \(\alpha\) around the vertical axis so that;

2. rotate through angle \(\beta\) around the lateral axis so that.

Mathematically this is given by:

()\[\begin{split}\begin{aligned} \left[ \begin{array}{c}{u_{2}} \\ {v_{2}} \\ {w_{2}}\end{array} \right]= \left[ \begin{array}{ccc}{\cos \alpha \cos \beta} & {\sin \alpha \cos \beta} & {\sin \beta} \\ {-\sin \alpha} & {\cos \alpha} & {0} \\ {-\cos \alpha \sin \beta} & {-\sin \alpha \sin \beta} & {\cos \beta}\end{array}\right] \times \left[\begin{array}{c} {\overline{u}} \\ {\overline{v}} \\ {\overline{w}} \end{array} \right] \end{aligned}\end{split}\]

where \(\alpha=\tan ^{-1}(\overline{v} / \overline{u})\) and \(\beta=\tan ^{-1}(\overline{w} / \sqrt{\overline{u}^{2}+\overline{v}^{2}})\).

Co-ordinate transformation of wind components

It is not generally possible to mount a 3-directional anemometer so that its axes coincide with the directions \(\overline{u}>0, \overline{v}=0, \overline{w}=0\) (where the over-bar denotes time-averaging over many data values). However, a co-ordinate transformation applied to the sensed components \(U, V, W\) means that the transformed component series \(u, v, w\) satisfies the above properties. The transformed components can be calculated using

()\[\begin{split}\begin{aligned} \begin{array}{l} {u=U \cos A \cos B+V \cos A \sin B+w \sin A} \\ {v=V \cos B-U \sin B} \\ {w=W \cos A-U \sin A \cos B-V \sin A \sin B} \end{array} \end{aligned}\end{split}\]

where angle \(A=\arctan{(\overline{W}/S)}\), \(\sin B=\overline{V}/S\), \(\sin B=\overline{U}/S\) , \(S=\overline{U}^2+\overline{V}^2\).

Errors in statistics

Reynolds averaging requires separation into a mean part (low frequency variation) and fluctuations (high frequency) from which we calculate covariances, variances, etc which are all in some sense a mean value taken over many samples. Usually we calculate the standard error of a mean.

()\[\Delta_{\alpha} = \sigma_{\alpha}/\sqrt{N}\]

where \(N\) is the number of samples taken. In turbulence, we know that each discrete measurement – or sample – is not fully independent of the last one, and the number of samples which are correlated is determined by the integral time-scale LT. So N should be replaced by the number of independent samples:

()\[N_{i} = T/L_{T}\]

where \(T\) is the period over which data is being averaged.

Statistical stationarity of a time series means that variances and covariances approach a stable value as the averaging time is extended, and the errors associated get smaller. So how long is long enough? The aim is to have a large number of samples.

So given that the averaging period \(T = N\Delta t\), there is a trade-off between the sampling period and the interval between samples, \(\Delta t\). If \(\Delta t\) is too long, then \(T\) must be increased to keep \(N\) large. The danger is that T is too long, and the statistics are no longer stationary, i.e. the turbulent flow has changed in response to external factors like a gust front passing through. Typically, sampling rate is 10-20 Hz, and the averaging period is 30-60 minutes depending on conditions.

The autocorrelation function and integral timescale

As well as calculating the covariance between two variables, it is instructive to look at the auto-correlation function, or the correlation of a variable with itself at later time-steps. For instance, consider the u component of the wind

()\[R_{u}(\tau) = \frac{\overline{u^{'}(t)u^{'}(t + \tau)}}{{\sigma_{u}}^{2}}\]

Hence \(R(\tau)=1\) at \(\tau = 0\). The rate at which \(R(\tau)\) decreases with lag is related to the size distribution of eddies. Large eddies cause slower variations in the time series, and thus the auto-correlation will decrease more slowly with lag than for a time series dominated by smaller eddies. Hence, a simple measure of ‘typical eddy size’ is given by the integral timescale \(L_T\) , defined as

()\[L_{T} = \int_{0}^{\infty}{R(\tau)\text{dτ}}\]

From Taylor’s frozen turbulence hypothesis, the integral length-scale \(L_X =\overline{u}L_T\). The integral length-scale for a variable can be thought of as the decorrelation length-scale, i.e. for two sensors separated beyond this distance, the turbulence measured at each will not be correlated.

Calculation of errors on covariances

For a covariance (e.g. \(\overline{w'T'})\ \) the error can be calculated: As \(\overline{w'T'}\) is a mean of a large number of quantities, we might expect its standard error to be given by:

()\[\Delta\overline{T^{'}w^{'}} = \sigma(T^{'}w^{'})/\sqrt{N_{i}}\]

where \(\sigma(T^{'}w^{'})\) is the standard deviation of \(\sigma(T^{'}w^{'})\) and Ni is the number of independent samples. Assuming all samples are independent of each other, estimate the standard error \(\Delta\overline{T^{'}w^{'}}\). Is this likely to be an accurate estimate? The number of independent samples is more accurately given by:

()\[N_{i} = N\frac{\Delta t}{L_{t}},\]

where \(N\) is the total number of samples, \(\Delta t\) is the time between samples and \(L_t\), is the integral timescale for the time series of \(T^{'}w^{'}\). To estimate the integral timescale, we need to calculate and plot the autocorrelation function for \(T^{'}w^{'}\), and read off the time at which it falls to \(1/e = 0.368\).

Note

Please report issues or ask questions about this site on the GitHub page.

Wind Profile

Wind profile in neutral conditions

The variation of mean wind speed \(\bar{u}\) with height z (in surface layer, above the RSL) above an aerodynamically rough surface can be represented by a logarithmic relation:

()\[\bar{u}(z) = \frac{u_{*}}{k}\ln\left\lbrack \frac{(z - d)}{z_{0}} \right\rbrack\]

where \(u_{*}\) is the friction velocity (\(u_{*}^{2} = - \overline{u'w'}\)) rate of vertical transfer by turbulence of horizontal momentum per unit mass of air), \(z_{0}\) is the roughness length of the surface for momentum, d is the zero-plane displacement and k is von Karman’s constant (0.4). The logarithmic law is strictly valid only in neutral conditions, i.e. when the effect of buoyancy on turbulence is small compared to the effect of wind shear. In such conditions, the temperature profile in the surface layer will be close to adiabatic (i.e. \(dT/dz=–9.8 \textrm{ K km}^{-1}\) or potential temperature is constant with height). When the sensible heat flux is significantly different from zero, Monin-Obukhov theory (or other) must be used.

Wind and temperature profile in non-neutral conditions

Modifications to the logarithmic profile are required in conditions of non-neutral stability. Here we use Monin-Obukhov theory of the surface layer that derives relations between the vertical variation of wind speed u(z) and potential temperature \(\theta(z)\) (which approximates the measured temperature T close to the surface), the scaling factors for momentum and temperature, \(u^*\) and \(T^*\), and the Monin‑Obukhov stability parameter

()\[\frac{z'}{L} = - \frac{k\left( z - d \right)\frac{g}{\theta_{0}}\frac{H}{\rho c_{p}}}{u_{*}^{3}}\]

where \(L\) is the Obukhov length and :math: \(z^’= z - d\). NB: the surface temperature \(\theta_0\) is an absolute temperature (units: K). The logarithmic profile relation can be rewritten for wind speed to include the stability corrections

()\[\bar{u}(z) = \frac{u_{*}}{k}\left\lbrack \ln\left( \frac{z - d}{z_{0}} \right) - \Psi_{m}\left( \frac{z - d}{L} \right) + \Psi_{m}\left( \frac{z_{0}}{L} \right) \right\rbrack\]

and similarly, for potential temperature:

()\[\bar{\theta}(z) = \theta_{0} + \frac{T_{*}}{k}\left\lbrack \ln\left( \frac{z - d}{z_{h}} \right) - \Psi_{h}\left( \frac{z - d}{L} \right) + \Psi_{h}\left( \frac{z_{h}}{L} \right) \right\rbrack\]

where the turbulent temperature scale is given by \(T_{*} = - \overline{w^{'}T^{'}}/u_{*} = - Q_{H}/(\rho c_{p}u_{*})\), \(\Psi_{m}\) is the integral stability correction function for momentum and \(\Psi_{h}\) is the integral stability correction function for heat. Note that both \(T_*\) and \(z’/L\) have the opposite sign to \(Q_H\) (which is positive in unstable conditions and negative in stable conditions). If \(z'/z \gg 1\) then the third term can assumed to be negligible (Garratt 1992)

Profiles and fluxes of moisture

Just as surface layer profiles of momentum and temperature follow a logarithmic form in neutral conditions, humidity in the surface layer has the same form, being transported by the same eddies. Thus, the profile is given by

()\[\overline{q}\left( z \right) - {\overline{q}}_{0} = \frac{q_{*}}{k}\ln\left( \frac{z - d}{z_{q}} \right)\]

where q is specific humidity, subscript 0 denotes a surface measurement, zq is the equivalent “roughness length” for humidity, and d is the zero-plane displacement height of a plant canopy over which measurements are being made). q* is a scaling variable, defined as

()\[q_{*} = \frac{- \overline{w'q'}}{u_{*}}\]

and thus the dimensionless moisture profile is defined as

()\[\phi_{q} = \frac{k\left( z - d \right)}{q_{*}}\frac{\partial\overline{q}}{\partial z}.\]

The moisture flux can be written in various equivalent forms

()\[E = \rho\overline{w'q'} = u_{*}q_{*} = - \rho K_{q}\frac{\partial\overline{q}}{\partial z}\]

where \(K_{q}\) is the eddy diffusivity for moisture. In neutral conditions it is assumed \(K_m =K_h =K_q=k(z-d)u_*\). Moisture follows Monin-Obukhov similarity just as other scalar variables do. This was not established at the Kansas experiments due to limitations in the accuracy of the measurements.

Note

Please report issues or ask questions about this site on the GitHub page.

Storage Heat Flux

The storage heat flux \(\Delta Q_S\) is the net energy stored or released by changes in sensible heat within the canopy air layer, roughness elements (e.g. vegetation, buildings in an urban environment), and the ground. Knowledge of \(\Delta Q_S\) is crucial to a wide range of processes and applications: from modelling turbulent heat transfer and boundary layer development to predicting soil thermal fields.

In rural sites, or simple bare soil sites, the flux may be a small fraction of the net all-wave radiation (Oliphant et al.,2004). However, in areas where there is more mass, such as cities, the term becomes much more significant (Kotthaus and Grimmond, 2014a) and a key element of the SEB and well-known effects such as the urban heat island.

Objective Hysteresis Model

The objective hysteresis model (OHM) is a simple model to calculate storage heat flux using net all-wave radiation \(Q^*\) following the hysteresis relation (Camuffo and Bernardi, 1982).

The OHM is forced by \(Q^∗\) and accounts for the diversity of the surface materials (sub-facets) in the measurement source area of interest with weightings (\(f\)) for their two- or three-dimensional extent (Grimmond et al., 1991):

()\[\Delta Q_{\mathrm{S}}=\sum_{i} f_{i}\left(a_{1, i} Q_{i}^{*}+a_{2, i} \frac{\partial Q_{i}^{*}}{\partial t}+a_{3, i}\right)\]

where the \(a_1\), \(a_2\), and \(a_3\) coefficients are for individual facets determined by least-square regression between \(\Delta Q_S\) and \(Q^*\) using results from observations.

These coefficients capture the net behaviour of a facet type in a typical setting, rather than being required to identify the component materials within a facet (e.g. multiple materials making up a roof, wall, with varying thermal connectivity and individual properties). A set of example OHM coefficients can be found here.

Note

Please report issues or ask questions about this site on the GitHub page.

Metrics for Model evaluation

Methods commonly used to evaluate model performance, include:

• Mean absolute error (MAE)

\[\mathrm{MAE}=\frac{1}{N} \sum_{i=1}^{N}\left|y_{i}-\hat{y}_{i}\right|\]

where \(N\) is number of observations, \(y_i\) the actual expected output and \(\hat{y}_{i}\) the model’s prediction (same notations below if not indicated otherwise).

• Mean bias error (MBE)

\[\mathrm{MBE}=\frac{1}{N} \sum_{i=1}^{N}\left(y_{i}-\hat{y}_{i}\right)\]

• Mean square error (MSE)

\[\mathrm{MSE}=\frac{1}{N} \sum_{i=1}^{N}\left(y_{i}-\hat{y}_{i}\right)^{2}\]

• Root mean square error (RMSE)

\[\mathrm{RMSE}=\sqrt{\frac{1}{N} \sum_{i=1}^{N}\left(y_{i}-\hat{y}_{i}\right)^{2}}\]

• Coefficient of determination (\(R^2\))

\[ \begin{align}\begin{aligned}R^{2}= 1-\frac{\mathrm{MSE}(\text { model })} {\mathrm{MSE}(\text { baseline })}\\\mathrm{MSE}(\text { baseline })= \frac{1}{N} \sum_{i=1}^{N}\left(y_{i}-\overline{y}\right)^{2}\end{aligned}\end{align} \]

where \(\overline{y}\) is mean of observed \(y_i\).

These presented with plots (e.g. scatter, time series) allow identification of periods where model perform well/poorly relative to observations. It should be remembered that both the model (e.g. parameters, forcing data) and the evaluation observations have errors.

Note

Please report issues or ask questions about this site on the GitHub page.

URAO

Safety

The Observatory should be regarded as a laboratory environment. The same guidelines apply as for other laboratory courses in the Department (see Departmental Handbook).

Note

1. Do not make adjustments to mains-supplied instruments or attempt to modify any apparatus without consulting either a demonstrator or a technician.

2. You should not work in the Observatory outside the class contact period, except by special arrangement with a supervising member of staff or a laboratory technician.

3. On the Observatory, beware of guy ropes attached to masts and overhead instruments. Be especially careful if the ground is muddy. Work in the area you have been shown to work in, as some cables and soil heat flux plates may be near the surface.

4. Wear suitable, protective footwear (i.e. NOT flip-flops) and beware of low-lying objects which you could trip over. Wear appropriately warm clothes.

5. Always ask if there are any points which are not clear.

Instruments

Overview

• Wind

Wind-related measurements at URAO

Instrument	Model	Manufacture	Sampling rate	comment
Ultrasonic anemometer	Windmaster pro	Gill	10 Hz	3 m
propeller anemometer	Gill	Gill	1 Hz	3 m
Wind profile anemometers	A100L2	Vector instruments	1 Hz	2,5, 10 m
Wind vane	W200P	Vector instruments	1 Hz	10 m
Vector Anemometer	A101M	Vector instruments	1 Hz	0.56, 0.8, 1.12, 1.6, 2.24, 3.2, 4.48, 6.4 m

• Heat and humidity

Temperature&humidity-related measurements at URAO

Instrument	Model	Manufacture	Sampling rate	comment
temperture probe	DTS12A	Vaisala	1 Hz	in Stevenson screen
Soil temperature	DTS12G	Vaisala	1 Hz	in soil (5cm-10cm-20cm-30cm-50cm-100cm)
humidity probe	HMP45A humidity probe	Vaisala	1 Hz	in Stevenson screen
platinum resistance themometer	PT100/3	Campbell scientific	1 Hz	temp profile (0.56m-1.12m-2.24m-4.48m)
Soil heat flux plate	HEP01	Hukseflux	1 Hz

• Radiation

Radiation-related measurements at URAO

Instrument	Model	Manufacture	Sampling rate	comment
pyranometer	CM11	Kipp and zonen	1 Hz	1m
net-radiometer	CNR4	Kipp and zonen	1 Hz	1.3m
Sunshine duration	CSD1	Kipp and zonen	1 Hz	10m

• Other

Other measurements at RUAO

Instrument	Model	Manufacture	Sampling rate	comment
infra-red gas analyser	Li-7500RS	Li-Cor	10 Hz	3m
Rain gauge	SBS 500	Environmental measurements Ltd.	1 Hz
Pressure	DPI1 140	Druck	1 Hz
Ceilometer	CL31	Vaisala	3sec row data	processed by R.Brugge 5min data

Eddy Covariance (EC) mast

The EC mast at the URAO has a sonic anemometer (Gill R3, Fig. 10.3) which derives the wind-speed from transit times of acoustic pulses travelling in both directions along a fixed path. The wind-speed component along the path is proportional to the difference between the transit times. Three sets of transducer pairs are used to derive the three components of the wind vector u, v, w. In addition, the speed of sound can be deduced:

Also, on the EC mast is an open path infra-red gas analyser (Li-Cor, LI-7500, Fig. 10.2) to measures both CO2 and H2O. It uses the principle of radiation absorption by water and carbon dioxide molecules at certain wavelengths. A differential measurement is made to compare transmittance at a wavelength with strong absorption adjacent to a wavelength where absorption is negligible. For water vapour (carbon dioxide) the wavelength is 2.59 \(\mathrm{\mu m}\) (4.26 \(\mathrm{\mu m}\)) and the reference wavelength is 2.4 \(\mathrm{\mu m}\) (3.95 \(\mathrm{\mu m}\)). The ratio of light intensity at the two wavelengths is proportional to the amount of water vapour present.

The IRGA at URAO is an open-path rather than a closed-path instrument (where air is sucked down a tube into the instrument itself). The q specific humidity of water vapour is expressed in units of kg kg^-1. The absolute humidity (\(\text{kg m}^{-3}\)) is derived by taking the molecular weight of water into account (1 mol = 18 g = 0.018 kg) and similarly for carbon dioxide concentrations (molar mass 44 g \(\text{mol}^{-1}\)). The instruments are mounted close to each other at a height of 3 m. A Campbell CR3000 logger is used to record the data at a sampling rate of 10 Hz.

Wind and temperature profile mast (6.4 m)

A profile of 8 pulse cup anemometers and 4 platinum resistance thermometers (PRTs) are mounted at various heights (Table 14.5). The coincident temperature and wind profiles allow both stability and surface fluxes to be derived. Each anemometer produces electrical pulses at a rate proportional to its rotation speed. The PRT output voltage is proportional to the PRT resistance.

measurement heights of temperature (T) and wind speed (U) at URAO

Variables measured	U&T	U	U&T	U	U&T	U	U&T	U
Height z (m)	0.56	0.80	1.12	1.60	2.24	3.20	4.48	6.40

Logging of sensors

Programmed data loggers sample the data at different time intervals. Raw samples (e.g. from EC system) or just statistics (e.g. an average from pyranometer) are recorded. During data processing calibration coefficients are applied.

Data from the Observatory

Data can be downloaded from: https://metdata.reading.ac.uk/ext/

Ask your instructor for download token in class if you need one.

Types of data

1. 5 min averaged logger output.

   • Includes individual radiation fluxes, soil heat flux, temperature (T), wind speed (WS), wind direction (Wdir), station pressure, rainfall, and relative humidity (RH).

2. Eddy covariances - 30-min averages.

   • Fully processed EC fluxes: These have been subjected to the numerous corrections (Kotthaus and Grimmond 2012, 2014a) that are regularly undertaken for EC fluxes.

3. 5-min WMO-standard processed output:

   • This includes the wind profile data and the temperature profile data. Radiation data (make certain you use the corrected longwave radiation data)

4. 0.1s Sonic Licor

   • Raw EC data - these files are very large so do NOT download data until you know what you really want/need.

References: See References.

Note

Please report issues or ask questions about this site on the GitHub page.

References

1

Arya, Pal S. Introduction to Micrometeorology. Volume 79. Elsevier, 2001. ISBN 9780120593545. URL: https://doi.org/10.1016/s0074-6142(01)x8014-5, doi:10.1016/s0074-6142(01)x8014-5.

2

Arya, Pal S. Air Pollution Meteorology. Elsevier, 2002. ISBN 9781898563938. URL: https://doi.org/10.1016/c2013-0-18095-6, doi:10.1016/c2013-0-18095-6.

3

Aubinet, Marc, Vesala, Timo, and Papale, Dario. Eddy Covariance. Springer Netherlands, 2012. ISBN 9789400723504, 9789400723511. URL: https://doi.org/10.1007/978-94-007-2351-1, doi:10.1007/978-94-007-2351-1.

4

Bradley, Stuart. Atmospheric acoustic remote sensing: Principles and applications. CRC press, 2007.

5

Burba, G. Eddy covariance method for scientific, industrial, agricultural, and regulatory applications : A field book on measuring ecosystem gas exchange and areal emission rates. LI-COR Biosciences, 2013. ISBN 061576827X 9780615768274. URL: https://bit.ly/35aOVD8.

6

Emeis, Stefan. Surface-Based Remote Sensing of the Atmospheric Boundary Layer. Volume 40. Springer Netherlands, 2011. ISBN 9789048193394, 9789048193400. URL: https://doi.org/10.1007/978-90-481-9340-0, doi:10.1007/978-90-481-9340-0.

7

Foken, Thomas. Micrometeorology. Volume 2. Springer Berlin Heidelberg, 2017. ISBN 9783642254390, 9783642254406. URL: https://doi.org/10.1007/978-3-642-25440-6, doi:10.1007/978-3-642-25440-6.

8

Garratt, J. R. The atmospheric boundary layer. Cambridge University Press, 1994. ISBN 9780521467452.

9

Grimmond, C.S.B. Aerodynamic roughness of urban areas derived from wind observations. Boundary Layer Meteorol., 89(1):1–24, October 1998. URL: https://doi.org/10.1023/a:1001525622213, doi:10.1023/a:1001525622213.

10

Grimmond, C.S.B., Cleugh, H.A., and Oke, T.R. An objective urban heat storage model and its comparison with other schemes. Atmospheric Environment. Part B. Urban Atmosphere, 25(3):311–326, January 1991. URL: https://doi.org/10.1016/0957-1272(91)90003-w, doi:10.1016/0957-1272(91)90003-w.

11

Kaimal, J. C. and Finnigan, J. J. Atmospheric Boundary Layer Flows : Their Structure and Measurement. Oxford University Press, Incorporated, 1994. ISBN 9780195362770.

12

Kotthaus, Simone and Grimmond, C.S.B. Identification of micro-scale anthropogenic CO2, heat and moisture sources – processing eddy covariance fluxes for a dense urban environment. Atmos. Environ., 57:301–316, September 2012. URL: https://doi.org/10.1016/j.atmosenv.2012.04.024, doi:10.1016/j.atmosenv.2012.04.024.

13

Kotthaus, Simone and Grimmond, C.S.B. Energy exchange in a dense urban environment – part i: Temporal variability of long-term observations in central London. Urban Clim., 10:261–280, December 2014. URL: https://doi.org/10.1016/j.uclim.2013.10.002, doi:10.1016/j.uclim.2013.10.002.

14

Leclerc, Monique Y. and Foken, Thomas. Footprints in Micrometeorology and Ecology. Volume 239. Springer Berlin Heidelberg, 2014. ISBN 9783642545443, 9783642545450. URL: https://doi.org/10.1007/978-3-642-54545-0, doi:10.1007/978-3-642-54545-0.

15

LeMone, Margaret A., Angevine, Wayne M., Bretherton, Christopher S., Chen, Fei, Dudhia, Jimy, Fedorovich, Evgeni, Katsaros, Kristina B., Lenschow, Donald H., Mahrt, Larry, Patton, Edward G., Sun, Jielun, Tjernström, Michael, and Weil, Jeffrey. 100 years of progress in boundary layer meteorology. Meteor. Monogr., 59():9.1–9.85, 2018. URL: https://doi.org/10.1175/amsmonographs-d-18-0013.1, doi:10.1175/amsmonographs-d-18-0013.1.

16

Moene, Arnold F. and van Dam, Jos C. Transport in the Atmosphere-Vegetation-Soil Continuum. Cambridge University Press, 2013. ISBN 9781139043137. URL: https://doi.org/10.1017/cbo9781139043137, doi:10.1017/cbo9781139043137.

17

Monson, Russell and Baldocchi, Dennis. Terrestrial Biosphere-Atmosphere Fluxes. Cambridge University Press, 2014. ISBN 9781139629218. URL: https://doi.org/10.1017/cbo9781139629218, doi:10.1017/cbo9781139629218.

18

Oke, T. R. Boundary Layer Climates. Routledge, September 2002. ISBN 9780203407219. URL: https://doi.org/10.4324/9780203407219, doi:10.4324/9780203407219.

19

Oke, T. R., Mills, G., Christen, A., and Voogt, J. A. Urban Climates. Cambridge University Press, 2017. ISBN 9781139016476. URL: https://doi.org/10.1017/9781139016476, doi:10.1017/9781139016476.

20

Stull, Roland B. An Introduction to Boundary Layer Meteorology. Volume 13. Springer Netherlands, 1988. ISBN 9789027727695, 9789400930278. URL: https://doi.org/10.1007/978-94-009-3027-8, doi:10.1007/978-94-009-3027-8.

21

Ward, H.C., Kotthaus, S., Järvi, L., and Grimmond, C.S.B. Surface urban energy and water balance scheme (SUEWS): Development and evaluation at two UK sites. Urban Clim., 18:1–32, December 2016. URL: https://doi.org/10.1016/j.uclim.2016.05.001, doi:10.1016/j.uclim.2016.05.001.

22

Wilczak, J. M., Gossard, E. E., Neff, W. D., and Eberhard, W. L. Ground-based remote sensing of the atmospheric boundary layer: 25 years of progress. In Garratt, J. R. and Taylor, P. A., editors, Boundary-Layer Meteorology 25th Anniversary Volume, 1970–1995, pages 321–349. Springer Netherlands, Dordrecht, 1996. URL: https://doi.org/10.1007/978-94-017-0944-6_14, doi:10.1007/978-94-017-0944-6_14.

23

Wyngaard, John C. Turbulence in the Atmosphere. Cambridge University Press, 2009. ISBN 9780511840524. URL: https://doi.org/10.1017/cbo9780511840524, doi:10.1017/cbo9780511840524.

Note

Please report issues or ask questions about this site on the GitHub page.

Python

Useful Links

There are a lot of excellent tutorials and web pages for helping with Python coding:

• The gist of Python: a quick introductory blog that covers Python basics for data analysis.

• Jupyter Notebook: Jupyter Notebook provides a powerful notebook-based data analysis environment that SuPy users are strongly encouraged to use. Jupyter notebooks can run in browsers (desktop, mobile) either by easy local configuration or on remote servers with pre-set environments (e.g., Google Colaboratory, Microsoft Azure Notebooks). In addition, Jupyter notebooks allow great shareability by incorporating source code and detailed notes in one place, which helps users to organise their computation work.

  • Installation

    Jupyter notebooks can be installed with pip on any desktop/server system and open .ipynb notebook files locally:

    python3 -m pip install jupyter -U
    

  • Extensions: To empower your Jupyter Notebook environment with better productivity, please check out the Unofficial Jupyter Notebook Extensions. Quick introductory blogs can be found here and here.

• pandas: pandas is an essential tool for data analysis in Python.

  • Introductory blogs:

    • Quick dive into Pandas for Data Science: introduction to pandas.

    • Basic Time Series Manipulation with Pandas: pandas-based time series manipulation.

    • Introduction to Data Visualization in Python: plotting using pandas and related libraries.

  • A detailed tutorial in Jupyter Notebooks:

    • Introduction to pandas

    • pandas fundamentals

    • Data Wrangling with pandas

Note

Please report issues or ask questions about this site on the GitHub page.

Jupyter Notebook

The Jupyter Notebook is an ideal tool for exploratory data analysis, by using which you can finish data processing, computation, plotting and even writing (using Markdown) in one stop. Also, it supports more programming languages than just Python: Julia, R, etc.. In this course, we will only use Python for its well developed ecosystem in data science and scientific computing.

Below we only cover several key points for using Jupyter Notebooks A more complete tutorial can be found here.

Jupyter Notebook Extensions

In addition to the Jupyter Notebook itself, we recommend Jupyter Notebook Extensions as well, which can further enhance your productivity in the Jupyter Notebook environment. A good guide can be found here.

Among a variety of extensions, the following are considered essential and particularly useful for this course, and many other data-analysis-oriented ones:

1. Table of Contents: this extension helps building well-structured notebooks and allows easier navigation.

2. Collapsible Headings: this tool further enhances the power of the above one by allowing collasible sections, particularly useful for large notebooks.

3. ExecuteTime: this gadget reports the executation time of each cell so you can always have a measure of the code efficiency.

4. Code prettify: this tool quickly formats your Python code with a pretty and tidy layout.

Best Practices

A good guide can be found here.

And we emphasise that structuring your notebook is the key. This post helps you set up your Notebook with a good structure.

Note

Please report issues or ask questions about this site on the GitHub page.

GitHub

For those who are new to GitHub, here is a quick introduction by GitHub themselves:

“GitHub is a code hosting platform for version control and collaboration. It lets you and others work together on projects from anywhere.”

A good tutorial for GitHub can be found here. Also, we suggest students to go through more tutorials below to get familiar with git and GitHub:

1. Git Handbook: Learn about version control—in particular, Git, and how it works with GitHub.

2. GitHub-based workflow: This guide explains how and why GitHub flow works.

Note

Please report issues or ask questions about this site on the GitHub page.

Python Tutorials

The following tutorials are provided to demonstrate basic skills in data analysis and visualisation using Python.

Note

These tutorials are in Jupyter notebooks and can be executed and viewed online.

This page was generated from /home/docs/checkouts/readthedocs.org/user_builds/blm/checkouts/latest/docs/tutorials/tutorial-1-basic-plotting.ipynb. Interactive online version: Slideshow:

Basic plotting in Python

Here you will learn the basics of plotting in python including:

1. How to have a basic plot?

2. How to add labels, legends and title to plot?

3. How to change the color, fontsize, and other properties of the lines in the plot?

4. How to save plots to files with different formats?

Loading necessary packages

import numpy as np
import matplotlib.pyplot as plt

Let’s create some sample data for plots

x=np.arange(0,10,.1)
y=np.sin(x)

Let’s have a quick plot:

plt.plot(x,y)

[<matplotlib.lines.Line2D at 0x112fd9748>]

We want to add x and y labels and title to the plot:

plt.plot(x,y)
plt.xlabel('This is x label')
plt.ylabel('This is y label')
plt.title('This is the title')

Text(0.5, 1.0, 'This is the title')

Now let’s add more lines to the plot, and try to make legends

plt.plot(x,y,label='y=sin(x)')
plt.plot(x,2*y,label='x=2*sin(x)')
plt.xlabel('This is x label')
plt.ylabel('This is y label')
plt.title('This is the title')
plt.legend()

<matplotlib.legend.Legend at 0x1131e93c8>

how to have more complex labels such as mathematic formula:

For math formulas, you need to use: $your formula here$

let’s give it a try:

plt.plot(x,x*x,label='$y=x^2$')
plt.plot(x,x*x*x,label='$y=x^3$')
plt.xlabel('This is x label')
plt.ylabel('This is y label')
plt.title('This is the title')
plt.legend()

<matplotlib.legend.Legend at 0x11324d630>

As you can see above, the superscript is made using ^{your superscript} and subscript is made using _{your subscript}

$x^{-3}$=\(x^{-3}\)

$x_{data}$=\(x_{data}\)

How to change the line width, line style and line color:

plt.plot(x,x*x,label='$y=x^2$',linewidth=2,color='r',linestyle='--')
plt.plot(x,x*x*x,label='$y=x^3$',linewidth=4,color='b')
plt.xlabel('This is x label')
plt.ylabel('This is y label')
plt.title('This is the title')
plt.legend()

<matplotlib.legend.Legend at 0x1133139b0>

How to change the text fontsize:

plt.plot(x,x*x,label='$y=x^2$',linewidth=2,color='r',linestyle='--')
plt.plot(x,x*x*x,label='$y=x^3$',linewidth=4,color='b')
plt.xlabel('this is fontsize 14',fontsize=14)
plt.ylabel('this is fontsize 20',fontsize=20)
plt.title('This is the title')
plt.legend(fontsize=14)

<matplotlib.legend.Legend at 0x1134b1eb8>

how to change the figure size:

plt.figure(figsize=(20,5))
plt.plot(x,x*x,label='$y=x^2$',linewidth=2,color='r',linestyle='--')
plt.plot(x,x*x*x,label='$y=x^3$',linewidth=4,color='b')
plt.xlabel('this is fontsize 14',fontsize=14)
plt.ylabel('this is fontsize 20',fontsize=20)
plt.title('This is the title')
plt.legend(fontsize=14)

<matplotlib.legend.Legend at 0x11356fe10>

and finally, how to save your figure to a file:

plt.figure(figsize=(20,5))
plt.plot(x,x*x,label='$y=x^2$',linewidth=2,color='r',linestyle='--')
plt.plot(x,x*x*x,label='$y=x^3$',linewidth=4,color='b')
plt.xlabel('this is fontsize 14',fontsize=14)
plt.ylabel('this is fontsize 20',fontsize=20)
plt.title('This is the title')
plt.legend(fontsize=14)
plt.savefig('sample.png',dpi=100)

Note that you can use various formats such as jpg, pdf, etc.

Also you can control the quality of the saved file using dpi keyword. Be careful that using large numbers for dpi would make the process of saving figures slow and the created file very large.

This page was generated from /home/docs/checkouts/readthedocs.org/user_builds/blm/checkouts/latest/docs/tutorials/tutorial-2-advance-plotting.ipynb. Interactive online version: Slideshow:

Advanced plotting in Python

Here you will learn more advance skills for plotting in Python:

1. How to have multiple plots in one figure (subplots)?

2. How to handle different axis in one figure?

3. How to position legend?

4. How to change x and y ticks?

Loading necessary packages

import numpy as np
import matplotlib.pyplot as plt

We want to have a figure with two horizontal subplots:

fig,axs=plt.subplots(1,2,figsize=(15,3))

The above figures have two axis (left and right one). You can see this by printing the contents of axs:

print(axs)

[<matplotlib.axes._subplots.AxesSubplot object at 0x115cbdc18>
 <matplotlib.axes._subplots.AxesSubplot object at 0x115d85b38>]

so let’s create some data to plot:

x=np.arange(0,10,.1)
y1=np.sin(x)
y2=np.cos(x)

Let’s plot them on the figure:

fig,axs=plt.subplots(1,2,figsize=(15,3))

ax1=axs[0] #first axis (left one)
ax2=axs[1] #second axis (right one)

ax1.plot(x,y1,color='r',label='left axis')
ax1.legend()

ax2.plot(x,y2,color='b',label='right axis')
ax2.legend()

<matplotlib.legend.Legend at 0x115a6d198>

As you can see above, it is as easy as defining on which axis you like to plot, and then everything is similar to single plots (almost everything, you see later on why). Now Let’s have more subplots:

fig,axs=plt.subplots(2,2,figsize=(15,3))

ax11=axs[0][0] # Top left
ax12=axs[0][1] # Top right
ax21=axs[1][0] # Bottom left
ax22=axs[1][1] # Bottom right

ax11.plot(x,y1,color='r',label='Top left')
ax11.legend()
ax11.set_ylabel('ylabel of top left plot')

ax12.plot(x,y2,color='b',label='Top right')
ax12.legend()


ax21.plot(x,y1,color='r',linestyle='--',label='Bottom left')
ax21.legend()

ax22.plot(x,y2,color='b',linestyle='--',label='Bottom right')
ax22.legend()
ax22.set_xlabel('xlabel of bottom right plot')

Text(0.5, 0, 'xlabel of bottom right plot')

Important note:

As you can see, in contrast to single plots, when you are trying to set a property for the plot, you use set_{keyword}.

For example:

in single plot: plt.xlabel('your xlabel')

in subplots plot: ax.set_xlabel('your xlabel')

If you look at the plot above, you can see that the subplots are positioned very intensly. To fix this, we use the following command: plt.tight_layout()

fig,axs=plt.subplots(2,2,figsize=(15,3))
plt.tight_layout()

ax11=axs[0][0] # Top left
ax12=axs[0][1] # Top right
ax21=axs[1][0] # Bottom left
ax22=axs[1][1] # Bottom right

ax11.plot(x,y1,color='r',label='Top left')
ax11.legend()
ax11.set_ylabel('ylabel of top left plot')

ax12.plot(x,y2,color='b',label='Top right')
ax12.legend()


ax21.plot(x,y1,color='r',linestyle='--',label='Bottom left')
ax21.legend()

ax22.plot(x,y2,color='b',linestyle='--',label='Bottom right')
ax22.legend()
ax22.set_xlabel('xlabel of bottom right plot')

Text(0.5, 6.000000000000025, 'xlabel of bottom right plot')

To adjust the distance between the subplots, you should use: fig.subplots_adjust(hspace=)

For example:

fig,axs=plt.subplots(2,2,figsize=(15,3))
fig.subplots_adjust(hspace=2)

ax11=axs[0][0] # Top left
ax12=axs[0][1] # Top right
ax21=axs[1][0] # Bottom left
ax22=axs[1][1] # Bottom right

ax11.plot(x,y1,color='r',label='Top left')
ax11.legend()
ax11.set_ylabel('ylabel of top left plot')

ax12.plot(x,y2,color='b',label='Top right')
ax12.legend()


ax21.plot(x,y1,color='r',linestyle='--',label='Bottom left')
ax21.legend()

ax22.plot(x,y2,color='b',linestyle='--',label='Bottom right')
ax22.legend()
ax22.set_xlabel('xlabel of bottom right plot')

Text(0.5, 0, 'xlabel of bottom right plot')

You can change the position of legend using the keyword loc=. Possible values for this keyword:

best
upper right
upper left
lower left
lower right
right
center left
center right
lower center
upper center
center

fig,axs=plt.subplots(1,1,figsize=(15,3))
axs.plot(x,y1,color='r',label='Top left')
axs.legend(loc='lower left',fontsize=16)

<matplotlib.legend.Legend at 0x11a22a3c8>

You can have costume x or y ticks using the following method:

fig,axs=plt.subplots(1,1,figsize=(15,3))
axs.plot(x,y1,color='r')
axs.set_xticks([0,2,4,6,8,10,12])

axs.set_xticklabels(['zero','two','four','six','eight','ten','twelve'])

[Text(0, 0, 'zero'),
 Text(0, 0, 'two'),
 Text(0, 0, 'four'),
 Text(0, 0, 'six'),
 Text(0, 0, 'eight'),
 Text(0, 0, 'ten'),
 Text(0, 0, 'twelve')]

This page was generated from /home/docs/checkouts/readthedocs.org/user_builds/blm/checkouts/latest/docs/tutorials/tutorial-3-Plotting albedo values in one plot.ipynb. Interactive online version: Slideshow:

Plotting albedo values in one plot

Here you can learn how to plot the albedo of two different days in one plot

First, some packages:

# These packages are necessary later on. Load all the packages in one place for consistency
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from pathlib import Path
import datetime

Let’s load the data and calculate albedos for the all data

#The path of the directory where all AMF data are
path_dir = Path.cwd()/'data'/'1'
name_of_site = 'CA-Obs_clean.csv.gz'
path_data = path_dir/name_of_site
path_data.resolve()
df_data = pd.read_csv(path_data, index_col='time',parse_dates=['time'])
ser_alb=df_data['SWOUT']/df_data['SWIN']

Let’s plot the albedo of one day

date1 = '2001 10 27'
ser_alb[ser_alb.between(0,1)&(df_data['SWIN']>5)].loc[date1].plot(marker='o')

<matplotlib.axes._subplots.AxesSubplot at 0x114ed4748>

And for another day

date2 = '2001 11 21'
ser_alb[ser_alb.between(0,1)&(df_data['SWIN']>5)].loc[date2].plot(marker='o')

<matplotlib.axes._subplots.AxesSubplot at 0x113672dd8>

Let’s first try to put the albedo of both days in one plot using previous methods

date1 = '2001 10 21'
ser_alb[ser_alb.between(0,1)&(df_data['SWIN']>5)].loc[date1].plot(marker='o')
date2 = '2001 11 21'
ser_alb[ser_alb.between(0,1)&(df_data['SWIN']>5)].loc[date2].plot(marker='o')

<matplotlib.axes._subplots.AxesSubplot at 0x11394f438>

You see that since the days are different, we do not get a good plot.

In order to fix this, we need to find the time of the day each data point is associated with. For example, for date1, we can find it with using the time level of the index as Dataframe.index.time :

ser_alb[ser_alb.between(0,1)&(df_data['SWIN']>5)].loc[date1].index.time

array([datetime.time(8, 30), datetime.time(9, 0), datetime.time(9, 30),
       datetime.time(10, 0), datetime.time(10, 30), datetime.time(11, 0),
       datetime.time(11, 30), datetime.time(12, 0), datetime.time(12, 30),
       datetime.time(13, 0), datetime.time(13, 30), datetime.time(14, 0),
       datetime.time(14, 30), datetime.time(15, 0), datetime.time(15, 30),
       datetime.time(16, 0), datetime.time(16, 30), datetime.time(17, 0),
       datetime.time(17, 30)], dtype=object)

Now let’s get this for both dates

X1=ser_alb[ser_alb.between(0,1)&(df_data['SWIN']>5)].loc[date1].index.time
X2=ser_alb[ser_alb.between(0,1)&(df_data['SWIN']>5)].loc[date2].index.time
Y1=ser_alb[ser_alb.between(0,1)&(df_data['SWIN']>5)].loc[date1]
Y2=ser_alb[ser_alb.between(0,1)&(df_data['SWIN']>5)].loc[date2]

Let’s try to plot them then in one figure

plt.plot(X1,Y1,'r',label='date 1')
plt.plot(X2,Y2,'b',label='date 2')
plt.legend()
plt.ylabel('albedo')

Text(0, 0.5, 'albedo')

Much better! but still the xticks are a little bit messy. We can use the follwoing to fix this issue

plt.plot(X1,Y1,'r',label='date 1')
plt.plot(X2,Y2,'b',label='date 2')
plt.legend()
plt.ylabel('albedo')
plt.xticks([datetime.time(x, 0) for x in [8,12,16,20]])

([<matplotlib.axis.XTick at 0x11430ba58>,
  <matplotlib.axis.XTick at 0x114383198>,
  <matplotlib.axis.XTick at 0x11440a7f0>,
  <matplotlib.axis.XTick at 0x114415630>],
 <a list of 4 Text xticklabel objects>)

This page was generated from /home/docs/checkouts/readthedocs.org/user_builds/blm/checkouts/latest/docs/tutorials/tutorial-4-Functions in Python.ipynb. Interactive online version: Slideshow:

Functions in Python

In this tutorial, you will learn how to use functions in Python and why are they extremely useful for the BLM course

What you will learn here:

1. How to write simple functions in Python?

2. How to call functions?

3. How to pass variables and parameters to functions?

4. How to use return in functions?

5. How to integrate functions together?

These packages are necessary later on. Load all the packages in one place for consistency

import numpy as np
import pandas as pd

Let’s say you want a function that simply prints Hello World for you. You can define this function as following:

def your_printer():
    print('Hello World')

As you can see, it is very simple to define the function, and when you run the above cell, nothing happens because we have not called the function yet. This is the handy part of functions in Python (or any other programming languages). You define the function once, and then you can all it whenever you want, and it does the job for you. You can call the function simply:

your_printer()

Hello World

Congrats! you just wrote and called the first function in Python.

But the application of functions goes beyond just printing. You can pass variables into function and they do the job for you. For example, let’s look modify your_printer function to print what we pass to it:

def your_printer(your_word):
    print(your_word)

Now let’s call this function:

your_printer('This is your word!')

This is your word!

You can call a function multiple times, and this is where functions save you time and space in coding:

your_printer('first word')
your_printer('second word')
your_printer('third word')

first word
second word
third word

Sometimes, you would like to pass a variable to a function, and the function does some mathematical operations and return the results to you for future use. In this case, you use return in functions. For example

def sum_two_number(number1,number2):

    return number1+number2

The above function, simply gets two number and returns the sum of them to you. You can use it as follows

sum_numbers=sum_two_number(1,2)
print(sum_numbers)

3

Then you can use the result of the function to do other stuff. For example:

print(sum_two_number(sum_numbers,4))

7

Now let’s have a more complex example. Let’s say you have the following equations:

\(q_1=x^2+x\)

\(q_2=\sin(q_1)+\frac{1}{q_1+1}\)

\(q_3=q_2*q_1+x^3\)

\(q_4=q_1+q_2+q_3+x\)

Now given a value for \(x\), you like to find the value of \(q_4\). This is where functions come handy. Let’s implemet this:

############################## q1

def cal_q1(x):

    return x**2+x

############################## q2

def cal_q2(q1):

    return np.sin(q1)+1/(q1+1)

############################## q3

def cal_q3(q1,q2,x):

    return q2*q1+x**3

############################## q4

def cal_q4(x):

    q1=cal_q1(x)
    q2=cal_q2(q1)
    q3=cal_q3(q1,q2,x)

    return q1+q2+q3+x

now we can calculate values of \(q_4\) given any arbitrary value for \(x\)

cal_q4(12)

1758.5598148460792

cal_q4(1)

7.727892280477045

cal_q4(5)

130.3710196531213

This page was generated from /home/docs/checkouts/readthedocs.org/user_builds/blm/checkouts/latest/docs/tutorials/tutorial-5-Scatter_plot_datetime.ipynb. Interactive online version: Slideshow:

Getting subset of time series + Scatter plot

Here you can learn how to:

- plot time series in points (scatter plot)
- Get subset of time series data and plot them

First, some packages:

# These packages are necessary later on. Load all the packages in one place for consistency
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from pathlib import Path
import datetime

Let’s load the data:

#The path of the directory where all AMF data are
path_dir = Path.cwd()/'data'/'1'
name_of_site = 'CA-Obs_clean.csv.gz'
path_data = path_dir/name_of_site
path_data.resolve()
df_data = pd.read_csv(path_data, index_col='time',parse_dates=['time'])
df_data.head()

	WS	RH	TA	PA	WD	P	SWIN	LWIN	SWOUT	LWOUT	NETRAD	H	LE
time
1997-01-01 00:30:00	2.988	73.036	-24.570	93.942	102.84	NaN	-0.14	213.39	0.01	214.83	-1.59	NaN	NaN
1997-01-01 01:00:00	2.671	73.146	-24.562	93.887	96.09	NaN	-0.03	216.73	0.03	215.03	1.64	NaN	NaN
1997-01-01 01:30:00	2.303	73.151	-24.431	93.934	112.34	NaN	0.24	223.46	0.02	215.91	7.77	NaN	NaN
1997-01-01 02:00:00	2.789	73.093	-24.379	93.917	109.16	NaN	-0.04	218.32	0.00	215.72	2.56	NaN	NaN
1997-01-01 02:30:00	2.274	73.140	-24.284	93.977	115.74	NaN	0.14	217.89	0.01	215.99	2.02	NaN	NaN

Let’s have a plot of Net radiation:

df_data['NETRAD'].plot(figsize=(12,5))

<matplotlib.axes._subplots.AxesSubplot at 0x126df3400>

If you like to focus in particular part of the data, you use the following feature of pandas package (DataFrame.loc(index1:index2)):

df_sub=df_data.loc['2002 07 01':'2002 07 30']
df_sub.head()

	WS	RH	TA	PA	WD	P	SWIN	LWIN	SWOUT	LWOUT	NETRAD	H	LE
time
2002-07-01 00:00:00	3.508	51.637	14.268	93.494	266.48	0.0	0.13	286.54	-0.20	361.24	-74.37	-3.451	0.175
2002-07-01 00:30:00	3.183	53.221	13.799	93.509	265.86	0.0	0.07	282.77	-0.23	358.10	-75.02	NaN	NaN
2002-07-01 01:00:00	3.759	55.664	13.182	93.506	266.07	0.0	-0.11	281.14	-0.14	356.13	-74.96	-14.770	1.628
2002-07-01 01:30:00	4.518	52.553	13.217	93.500	276.19	0.0	0.02	282.81	0.15	362.15	-79.47	-69.300	12.360
2002-07-01 02:00:00	3.940	52.343	13.150	93.484	275.58	0.0	-0.35	282.55	0.15	362.48	-80.43	-35.270	3.713

df_sub.index

DatetimeIndex(['2002-07-01 00:00:00', '2002-07-01 00:30:00',
               '2002-07-01 01:00:00', '2002-07-01 01:30:00',
               '2002-07-01 02:00:00', '2002-07-01 02:30:00',
               '2002-07-01 03:00:00', '2002-07-01 03:30:00',
               '2002-07-01 04:00:00', '2002-07-01 04:30:00',
               ...
               '2002-07-30 19:00:00', '2002-07-30 19:30:00',
               '2002-07-30 20:00:00', '2002-07-30 20:30:00',
               '2002-07-30 21:00:00', '2002-07-30 21:30:00',
               '2002-07-30 22:00:00', '2002-07-30 22:30:00',
               '2002-07-30 23:00:00', '2002-07-30 23:30:00'],
              dtype='datetime64[ns]', name='time', length=1440, freq=None)

as you can see, we have a subset of df_data which contains the data from July 2002. Now let’s plot it:

df_sub['NETRAD'].plot(figsize=(15,5))
plt.ylabel('Temp ($\degree$C)')

Text(0, 0.5, 'Temp ($\\degree$C)')

You can even specify the date and time of interest in the subset:

df_sub.loc['2002 07 12 6:00:00':'2002 07 12 22:00:00']['TA'].plot(figsize=(15,5))
plt.ylabel('Temp ($\degree$C)')

Text(0, 0.5, 'Temp ($\\degree$C)')

Sometimes we want to plot points instead of lines, specially for time series when you have a lot of missing points. You can do this by using plt.scatter as follows:

Y=df_sub['TA']
X=df_sub.index

fig,ax=plt.subplots(1,1,figsize=(15,5))
plt.scatter(X,Y)
plt.xlim([df_sub.index.date.min(),df_sub.index.date.max()])
plt.ylabel('Temp ($\degree$C)')

Text(0, 0.5, 'Temp ($\\degree$C)')

You can change the marker color and style (For different types of markers, please refer here):

Y=df_sub['TA']
X=df_sub.index

fig,ax=plt.subplots(1,1,figsize=(15,5))
plt.scatter(X,Y,color='r',marker='+')
plt.xlim([df_sub.index.date.min(),df_sub.index.date.max()])
plt.ylabel('Temp ($\degree$C)')

Text(0, 0.5, 'Temp ($\\degree$C)')

This page was generated from /home/docs/checkouts/readthedocs.org/user_builds/blm/checkouts/latest/docs/tutorials/tutorial-AMF-sim.ipynb. Interactive online version: Slideshow:

Modelling Surface Energy Balance at an AmeriFlux Site Using SuPy

This tutorial aims to demonstrate how to use an advanced land surface model (SuPy, SUEWS in Python) to better understand the surface energy balance (SEB) features by conducting simulation at an AmeriFlux site. This would be particularly useful after building your own model: as you will learn how sophisticated models could be developed from those simpler ones.

SuPy is a Python-enhanced urban climate model with SUEWS, *Surface Urban Energy and Water Balance Scheme*, as its computation core. More SuPy tutorials are available here.

In this tutorial the workflow to model the surface energy balance (SEB) at a chosen AmeriFlux (AMF) site using SuPy/SUEWS is undertaken. The steps, consist of

1. Preparing the input data;

2. Running a simulation;

3. Examination of results; and

4. Further exploration

Before starting, you need to install SuPy and load the following necessary packages.

pip install supy==2019.11.18.dev0

# !pip install supy==2019.11.18.dev0 &> install.log

import matplotlib.pyplot as plt
import supy as sp
import pandas as pd
import numpy as np
from pathlib import Path
%matplotlib inline

%load_ext autoreload
%autoreload 2

sp.show_version()

supy: 2019.11.18dev
supy_driver: 2019a18

Prepare input data

Overview of SuPy input

Load sample data:

To ease the preparation of model input, a helper function load_SampleData is provided to get the sample input for SuPy simulations, which will later be used as template to populate your specific model configurations and forcing input.

df_state_init, df_forcing = sp.load_SampleData()

2019-11-22 09:16:28,573 — SuPy — INFO — All cache cleared.

df_state_init

df_state_init includes model Initial state consisting of:

• surface characteristics (e.g., albedo, emissivity, land cover fractions, etc.; full details refer to SUEWS documentation).

• model configurations (e.g., stability; full details refer to SUEWS documentation).

Detailed description of variables in df_state_init refers to SuPy input.

• Surface land cover fraction information in the sample input dataset:

df_state_init.filter(like='sfr')

var	sfr
ind_dim	(0,)	(1,)	(2,)	(3,)	(4,)	(5,)	(6,)
grid
98	0.43	0.38	0.001	0.019	0.029	0.001	0.14

• Heights of bluff-bodies (m):

df_state_init.loc[:, ['bldgh', 'evetreeh', 'dectreeh']]

var	bldgh	dectreeh	evetreeh
ind_dim	0	0	0
grid
98	22.0	13.1	13.1

df_forcing

df_forcing includes meteorological and other external forcing information.

Detailed description of variables in df_forcing refers to SuPy input.

Below is a view of heading lines of the forcing variables.

df_forcing.head()

	iy	id	it	imin	qn	qh	qe	qs	qf	U	...	snow	ldown	fcld	Wuh	xsmd	lai	kdiff	kdir	wdir	isec
2012-01-01 00:05:00	2012	1	0	5	-999.0	-999.0	-999.0	-999.0	-999.0	4.5225	...	-999.0	-999.0	-999.0	-999.0	-999.0	-999.0	-999.0	-999.0	-999.0	0.0
2012-01-01 00:10:00	2012	1	0	10	-999.0	-999.0	-999.0	-999.0	-999.0	4.5225	...	-999.0	-999.0	-999.0	-999.0	-999.0	-999.0	-999.0	-999.0	-999.0	0.0
2012-01-01 00:15:00	2012	1	0	15	-999.0	-999.0	-999.0	-999.0	-999.0	4.5225	...	-999.0	-999.0	-999.0	-999.0	-999.0	-999.0	-999.0	-999.0	-999.0	0.0
2012-01-01 00:20:00	2012	1	0	20	-999.0	-999.0	-999.0	-999.0	-999.0	4.5225	...	-999.0	-999.0	-999.0	-999.0	-999.0	-999.0	-999.0	-999.0	-999.0	0.0
2012-01-01 00:25:00	2012	1	0	25	-999.0	-999.0	-999.0	-999.0	-999.0	4.5225	...	-999.0	-999.0	-999.0	-999.0	-999.0	-999.0	-999.0	-999.0	-999.0	0.0

5 rows × 25 columns

Site-specific configuration of surface parameters

Given pandas.DataFrame as the core data structure of SuPy, all operations, including modification, output, demonstration, etc., on SuPy inputs (df_state_init and df_forcing) can be done using pandas-based functions/methods. Please see SuPy quickstart for methods to do so.

Below we will modify several key properties of the chosen site with appropriate values to run SuPy. First, we copy the df_state_init to have a new DataFrame for manipulation.

df_state_amf = df_state_init.copy()

# site identifier
name_site = 'US-AR1'

location

# latitude
df_state_amf.loc[:, 'lat'] = 41.37
# longitude
df_state_amf.loc[:, 'lng'] = -106.24
# altitude
df_state_amf.loc[:, 'alt'] = 611.

land cover fraction

Land covers in SUEWS

# view the surface fraction variable: `sfr`
df_state_amf.loc[:, 'sfr'] = 0
df_state_amf.loc[:, ('sfr', '(4,)')] = 1
df_state_amf.loc[:, 'sfr']

ind_dim	(0,)	(1,)	(2,)	(3,)	(4,)	(5,)	(6,)
grid
98	0.0	0.0	0.0	0.0	1.0	0.0	0.0

albedo

# we only set values for grass as the modelled site has a single land cover type: grass.
df_state_amf.albmax_grass = 0.19
df_state_amf.albmin_grass = 0.14

# initial albedo value
df_state_amf.loc[:, 'albgrass_id'] = 0.14

LAI/phenology

df_state_amf.filter(like='lai')

var	laimax			laimin			laipower								laitype			laicalcyes	lai_id
ind_dim	(0,)	(1,)	(2,)	(0,)	(1,)	(2,)	(0, 0)	(0, 1)	(0, 2)	(1, 0)	...	(3, 0)	(3, 1)	(3, 2)	(0,)	(1,)	(2,)	0	(0,)	(1,)	(2,)
grid
98	5.1	5.5	5.9	4.0	1.0	1.6	0.03	0.03	0.03	0.0005	...	0.0005	0.0005	0.0005	0.0	0.0	0.0	1	4.0	1.0	1.6

1 rows × 25 columns

# properties to control vegetation phenology
# you can skip the details for and just set them as provided below

# LAI paramters
df_state_amf.loc[:, ('laimax', '(2,)')] = 1
df_state_amf.loc[:, ('laimin', '(2,)')] = 0.2
# initial LAI
df_state_amf.loc[:, ('lai_id', '(2,)')] = 0.2

# BaseT
df_state_amf.loc[:, ('baset', '(2,)')] = 5
# BaseTe
df_state_amf.loc[:, ('basete', '(2,)')] = 20

# SDDFull
df_state_amf.loc[:, ('sddfull', '(2,)')] = -1000
# GDDFull
df_state_amf.loc[:, ('gddfull', '(2,)')] = 1000

surface resistance

# parameters to model surface resistance
df_state_amf.maxconductance = 18.7
df_state_amf.g1 = 1
df_state_amf.g2 = 104.215
df_state_amf.g3 = 0.424
df_state_amf.g4 = 0.814
df_state_amf.g5 = 36.945
df_state_amf.g6 = 0.025

measurement height

# height where forcing variables are measured/collected
df_state_amf.z = 2.84

urban feature

# disable anthropogenic heat by setting zero population
df_state_amf.popdensdaytime = 0
df_state_amf.popdensnighttime = 0

check df_state

# this procedure is to double-check proper values are set in `df_state_amf`
sp.check_state(df_state_amf)

2019-11-22 09:16:30,954 — SuPy — INFO — SuPy is validating `df_state`...
2019-11-22 09:16:31,120 — SuPy — INFO — All checks for `df_state` passed!

prepare forcing conditions

Here we use the a SuPy utility function read_forcing to read in forcing data from an external file in the format of SUEWS input. Also note, this read_forcing utility will also resample the forcing data to a proper temporal resolution to run SuPy/SUEWS, which is usually 5 min (300 s).

load and resample forcing data

# load forcing data from an external file and resample to a resolution of 300 s.
# Note this dataset has been gap-filled.
df_forcing_amf = sp.util.read_forcing('./data/US-AR1_2010_data_60.txt',
                                      tstep_mod=300)

# this procedure is to double-check proper forcing values are set in `df_forcing_amf`
_ = sp.check_forcing(df_forcing_amf)

2019-11-22 09:16:32,059 — SuPy — INFO — SuPy is validating `df_forcing`...
2019-11-22 09:16:34,374 — SuPy — ERROR — Issues found in `df_forcing`:
`kdown` should be between [0, 1400] but `-1.3057500000000002` is found at 2010-01-01 00:05:00

The checker detected invalid values in variable kdown: negative incoming solar radiation is found. We then need to fix this as follows:

# modify invalid values
df_forcing_amf.kdown = df_forcing_amf.kdown.where(df_forcing_amf.kdown > 0, 0)

# check `df_forcing` again
_ = sp.check_forcing(df_forcing_amf)

2019-11-22 09:16:34,466 — SuPy — INFO — SuPy is validating `df_forcing`...
2019-11-22 09:16:36,946 — SuPy — INFO — All checks for `df_forcing` passed!

examine forcing data

We can examine the forcing data:

list_var_forcing = [
    'kdown',
    'Tair',
    'RH',
    'pres',
    'U',
    'rain',
]
dict_var_label = {
    'kdown': 'Incoming Solar\n Radiation ($ \mathrm{W \ m^{-2}}$)',
    'Tair': 'Air Temperature ($^{\circ}}$C)',
    'RH': r'Relative Humidity (%)',
    'pres': 'Air Pressure (hPa)',
    'rain': 'Rainfall (mm)',
    'U': 'Wind Speed (m $\mathrm{s^{-1}}$)'
}
df_plot_forcing_x = df_forcing_amf.loc[:, list_var_forcing].copy().shift(
    -1).dropna(how='any')
df_plot_forcing = df_plot_forcing_x.resample('1h').mean()
df_plot_forcing['rain'] = df_plot_forcing_x['rain'].resample('1h').sum()

axes = df_plot_forcing.plot(
    subplots=True,
    figsize=(8, 12),
    legend=False,
)
fig = axes[0].figure
fig.tight_layout()
fig.autofmt_xdate(bottom=0.2, rotation=0, ha='center')
for ax, var in zip(axes, list_var_forcing):
    _ = ax.set_ylabel(dict_var_label[var])

Run simulations

Once met-forcing (via df_forcing_amf) and initial conditions (via df_state_amf) are loaded in, we call sp.run_supy to conduct a SUEWS simulation, which will return two pandas DataFrames: df_output and df_state_final.

df_output, df_state_final = sp.run_supy(df_forcing_amf, df_state_amf)

2019-11-22 09:16:40,809 — SuPy — INFO — ====================
2019-11-22 09:16:40,810 — SuPy — INFO — Simulation period:
2019-11-22 09:16:40,811 — SuPy — INFO —   Start: 2010-01-01 00:05:00
2019-11-22 09:16:40,811 — SuPy — INFO —   End: 2011-01-01 00:00:00
2019-11-22 09:16:40,812 — SuPy — INFO —
2019-11-22 09:16:40,813 — SuPy — INFO — No. of grids: 1
2019-11-22 09:16:40,813 — SuPy — INFO — SuPy is running in serial mode
2019-11-22 09:16:54,304 — SuPy — INFO — Execution time: 13.5 s
2019-11-22 09:16:54,305 — SuPy — INFO — ====================

df_output

df_output is an ensemble output collection of major SUEWS output groups, including:

• SUEWS: the essential SUEWS output variables

• DailyState: variables of daily state information

• snow: snow output variables (effective when snowuse = 1 set in df_state_init)

• RSL: profile of air temperature, humidity and wind speed within roughness sub-layer.

Detailed description of variables in df_output refers to SuPy output

df_output.columns.levels[0]

Index(['SUEWS', 'snow', 'RSL', 'DailyState'], dtype='object', name='group')

df_state_final

df_state_final is a DataFrame for holding:

1. all model states if save_state is set to True when calling sp.run_supy and supy may run significantly slower for a large simulation;

2. or, only the final state if save_state is set to False (the default setting) in which mode supy has a similar performance as the standalone compiled SUEWS executable.

Entries in df_state_final have the same data structure as df_state_init and can thus be used for other SUEWS simulations staring at the timestamp as in df_state_final.

Detailed description of variables in df_state_final refers to SuPy output

df_state_final.T.head()

	datetime	2010-01-01 00:05:00	2011-01-01 00:05:00
	grid	98	98
var	ind_dim
ah_min	(0,)	15.0	15.0
	(1,)	15.0	15.0
ah_slope_cooling	(0,)	2.7	2.7
	(1,)	2.7	2.7
ah_slope_heating	(0,)	2.7	2.7

Examine results

Thanks to the functionality inherited from pandas and other packages under the PyData stack, compared with the standard SUEWS simulation workflow, supy enables more convenient examination of SUEWS results by statistics calculation, resampling, plotting (and many more).

Ouptut structure

df_output is organised with MultiIndex (grid,timestamp) and (group,varaible) as index and columns, respectively.

df_output.head()

	group	SUEWS										...	DailyState
	var	Kdown	Kup	Ldown	Lup	Tsurf	QN	QF	QS	QH	QE	...	DensSnow_Paved	DensSnow_Bldgs	DensSnow_EveTr	DensSnow_DecTr	DensSnow_Grass	DensSnow_BSoil	DensSnow_Water	a1	a2	a3
grid	datetime
98	2010-01-01 00:05:00	0.0	0.0	265.638676	305.413842	-1.587667	-39.775166	0.0	-50.989269	11.054221	0.159883	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
	2010-01-01 00:10:00	0.0	0.0	265.638676	305.413842	-1.587667	-39.775166	0.0	-50.729902	10.795477	0.159259	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
	2010-01-01 00:15:00	0.0	0.0	265.638676	305.413842	-1.587667	-39.775166	0.0	-50.481342	10.547515	0.158661	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
	2010-01-01 00:20:00	0.0	0.0	265.638676	305.413842	-1.587667	-39.775166	0.0	-50.243138	10.309886	0.158086	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
	2010-01-01 00:25:00	0.0	0.0	265.638676	305.413842	-1.587667	-39.775166	0.0	-50.014860	10.082159	0.157534	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN

5 rows × 345 columns

Here we demonstrate several typical scenarios for SUEWS results examination.

The essential SUEWS output collection is extracted as a separate variable for easier processing in the following sections. More advanced slicing techniques are available in pandas documentation.

grid = df_state_amf.index[0]
df_output_suews = df_output.loc[grid, 'SUEWS']

Statistics Calculation

We can use .describe() method for a quick overview of the key surface energy balance budgets.

df_output_suews.loc[:, ['QN', 'QS', 'QH', 'QE', 'QF']].describe()

var	QN	QS	QH	QE	QF
count	105120.000000	105120.000000	105120.000000	105120.000000	105120.0
mean	118.291222	9.117453	53.611159	56.302547	0.0
std	213.150698	82.433076	72.044586	93.959856	0.0
min	-98.895138	-103.332153	-136.770097	-10.592811	0.0
25%	-32.252113	-46.922708	9.790903	0.633106	0.0
50%	-1.216564	-33.883630	25.756834	4.240341	0.0
75%	247.458581	55.231281	75.382465	67.207728	0.0
max	746.187700	262.700179	403.492712	445.496829	0.0

Plotting

Basic example

Plotting is very straightforward via the .plot method bounded with pandas.DataFrame. Note the usage of loc for to slices of the output DataFrame.

# a dict for better display variable names
dict_var_disp = {
    'QN': '$Q^*$',
    'QS': r'$\Delta Q_S$',
    'QE': '$Q_E$',
    'QH': '$Q_H$',
    'QF': '$Q_F$',
    'Kdown': r'$K_{\downarrow}$',
    'Kup': r'$K_{\uparrow}$',
    'Ldown': r'$L_{\downarrow}$',
    'Lup': r'$L_{\uparrow}$',
    'Rain': '$P$',
    'Irr': '$I$',
    'Evap': '$E$',
    'RO': '$R$',
    'TotCh': '$\Delta S$',
}

Peek at the simulation results:

grid = df_state_init.index[0]

ax_output = df_output_suews\
    .loc['2010-06-01':'2010-06-07',
         ['QN', 'QS', 'QE', 'QH', 'QF']]\
    .rename(columns=dict_var_disp)\
    .plot()
_ = ax_output.set_xlabel('Date')
_ = ax_output.set_ylabel('Flux ($ \mathrm{W \ m^{-2}}$)')
_ = ax_output.legend()

Plotting after resampling

The suggested runtime/simulation frequency of SUEWS is 300 s, which usually results a large output and may be over-weighted for storage and analysis. Also, you may feel apparent slowdown in producing the above figure as a large amount of data were used for the plotting. To slim down the result size for analysis and output, we can resample the default output very easily.

rsmp_1d = df_output_suews.resample('1d')
# daily mean values
df_1d_mean = rsmp_1d.mean()
# daily sum values
df_1d_sum = rsmp_1d.sum()

We can then re-examine the above energy balance at hourly scale and plotting will be significantly faster.

# energy balance
ax_output = df_1d_mean\
    .loc[:, ['QN', 'QS', 'QE', 'QH', 'QF']]\
    .rename(columns=dict_var_disp)\
    .plot(
            figsize=(10, 3),
            title='Surface Energy Balance',
        )
_ = ax_output.set_xlabel('Date')
_ = ax_output.set_ylabel('Flux ($ \mathrm{W \ m^{-2}}$)')
_ = ax_output.legend()

Then we use the hourly results for other analyses.

# radiation balance
ax_output = df_1d_mean\
    .loc[:, ['QN', 'Kdown', 'Kup', 'Ldown', 'Lup']]\
    .rename(columns=dict_var_disp)\
    .plot(
        figsize=(10, 3),
        title='Radiation Balance',
    )
_ = ax_output.set_xlabel('Date')
_ = ax_output.set_ylabel('Flux ($ \mathrm{W \ m^{-2}}$)')
_ = ax_output.legend()

# water balance
ax_output = df_1d_sum\
    .loc[:, ['Rain', 'Irr', 'Evap', 'RO', 'TotCh']]\
    .rename(columns=dict_var_disp)\
    .plot(
        figsize=(10, 3),
        title='Surface Water Balance',
    )
_ = ax_output.set_xlabel('Date')
_ = ax_output.set_ylabel('Water amount (mm)')
_ = ax_output.legend()

Get an overview of partitioning in energy and water balance at monthly scales:

# get a monthly Resampler
df_plot = df_output_suews.copy()
df_plot.index = df_plot.index.set_names('Month')
rsmp_1M = df_plot\
    .shift(-1)\
    .dropna(how='all')\
    .resample('1M', kind='period')
# mean values
df_1M_mean = rsmp_1M.mean()
# sum values
df_1M_sum = rsmp_1M.sum()

# month names
name_mon = [x.strftime('%b') for x in rsmp_1M.groups]
# create subplots showing two panels together
fig, axes = plt.subplots(2, 1, sharex=True)
# surface energy balance
_=df_1M_mean\
    .loc[:, ['QN', 'QS', 'QE', 'QH', 'QF']]\
    .rename(columns=dict_var_disp)\
    .plot(
        ax=axes[0],  # specify the axis for plotting
        figsize=(10, 6),  # specify figure size
        title='Surface Energy Balance',
        kind='bar',
    )
# surface water balance
_=df_1M_sum\
    .loc[:, ['Rain', 'Irr', 'Evap', 'RO', 'TotCh']]\
    .rename(columns=dict_var_disp)\
    .plot(
        ax=axes[1],  # specify the axis for plotting
        title='Surface Water Balance',
        kind='bar'
    )

# annotations
_ = axes[0].set_ylabel('Mean Flux ($ \mathrm{W \ m^{-2}}$)')
_ = axes[0].legend()
_ = axes[1].set_xlabel('Month')
_ = axes[1].set_ylabel('Total Water Amount (mm)')
_ = axes[1].xaxis.set_ticklabels(name_mon, rotation=0)
_ = axes[1].legend()

Save results to external files

The supy output can be saved as txt files for further analysis using supy function save_supy.

list_path_save = sp.save_supy(df_output, df_state_final)

for file_out in list_path_save:
    print(file_out.name)

98_2010_SUEWS_5.txt
98_2010_snow_5.txt
98_2010_RSL_5.txt
98_2010_DailyState.txt
98_2010_SUEWS_60.txt
98_2010_snow_60.txt
98_2010_RSL_60.txt
df_state.csv

More explorations into simulation results

In this section, we will use the simulation results to explore more features revealed by SuPy/SUEWS simulations but unavailable in your simple model.

Dynamics in rainfall and soil moisture deficit (SMD)

df_dailystate = df_output.loc[grid, 'DailyState'].dropna(
    how='all').resample('1d').mean()

# daily rainfall
ser_p = df_dailystate.P_day.rename('Rainfall')
ser_smd = df_output_suews.SMD
ser_smd_dmax = ser_smd.resample('1d').max().rename('SMD')

ax = pd.concat([ser_p, ser_smd_dmax], axis=1).plot(secondary_y='SMD',
                                                   figsize=(9, 4))
_ = ax.set_xlabel('Time (month)')

Variability in albedo

How does albedo change over time?

ser_alb = df_dailystate.AlbGrass
ax = ser_alb.plot()
_ = ax.set_xlabel('Time (month)')

How is albedo associated with vegetation phenology?

ser_lai = df_dailystate.LAI_Grass
pd.concat([ser_lai, ser_alb], axis=1).plot(secondary_y='AlbGrass',
                                           figsize=(9, 4))
ax = ser_lai.plot()
_ = ax.set_xlabel('Time (month)')

<matplotlib.axes._subplots.AxesSubplot at 0x7fc9989c9ac8>

ax_alb_lai = df_dailystate[['LAI_Grass', 'AlbGrass']].plot.scatter(
    x='LAI_Grass',
    y='AlbGrass',
)
ax_alb_lai.set_aspect('auto')

Variability in surface resistance

How does surface resistance vary over time?

ser_rs = df_output_suews.RS

• intra-annual

ax = ser_rs.resample('1d').median().plot()
_ = ax.set_xlabel('Time (month)')

• intra-daily

# a winter day
ax = ser_rs.loc['2010-01-22'].between_time('0830', '1600').plot()
_ = ax.set_xlabel('Time')

# a summer day
ax = ser_rs.loc['2010-07-01'].between_time('0530', '1900').plot()
_ = ax.set_xlabel('Time')

How is surface resistance associated with other surface properties?

# SMD
ser_smd = df_output_suews.SMD
df_x = pd.concat([ser_smd, ser_rs],
                 axis=1).between_time('1000', '1600').resample('1d').mean()
df_x = df_x.loc[df_x.RS < 500]
_ = df_x.plot.scatter(
    x='SMD',
    y='RS',
)

# LAI
df_x = pd.concat(
    [ser_lai,
     ser_rs.between_time('1000', '1600').resample('1d').mean()],
    axis=1)
df_x = df_x.loc[df_x.RS < 500]
_ = df_x.plot.scatter(
    x='LAI_Grass',
    y='RS',
)

How is surface resistance dependent on meteorological conditions?

cmap_sel = plt.cm.get_cmap('RdBu', 12)

# solar radiation
# colour by season
ser_kdown = df_forcing_amf.kdown
df_x = pd.concat([ser_kdown, ser_rs], axis=1).between_time('1000', '1600')
df_x = df_x.loc[df_x.RS < 1500]
df_plot = df_x.iloc[::20]
ax = df_plot.plot.scatter(x='kdown',
                          y='RS',
                          c=df_plot.index.month,
                          cmap=cmap_sel,
                          sharex=False)
fig = ax.figure
_ = fig.axes[1].set_title('month')
fig.tight_layout()

# air temperature
ser_ta = df_forcing_amf.Tair
df_x = pd.concat([ser_ta, ser_rs], axis=1).between_time('1000', '1600')
df_x = df_x.loc[df_x.RS < 1500]
df_plot = df_x.iloc[::15]
ax = df_plot.plot.scatter(x='Tair',
                          y='RS',
                          c=df_plot.index.month,
                          cmap=cmap_sel,
                          sharex=False)
fig = ax.figure
_ = fig.axes[1].set_title('month')
fig.tight_layout()

# air humidity
ser_rh = df_forcing_amf.RH
df_x = pd.concat([ser_rh, ser_rs], axis=1).between_time('1000', '1600')
df_x = df_x.loc[df_x.RS < 1500]
df_plot = df_x.iloc[::15]
ax = df_plot.plot.scatter(x='RH',
                          y='RS',
                          c=df_plot.index.month,
                          cmap=cmap_sel,
                          sharex=False)
fig = ax.figure
_ = fig.axes[1].set_title('month')
fig.tight_layout()

• Task:

Based on the above plots showing RS vs. met. conditions, explore these relationships again at the intra-daily scales.

Note

Please report issues or ask questions about this site on the GitHub page.

Activities

These activities are designed to be in a logical calculation order to help you understand different physical properties and processes.

Tip

The tasks designed here are based on Python 3. We recommend this software bundle (with download link):

1. Anaconda 3: an easy-to-use Python 3 bundle with many preinstalled data analysis packages. download Anaconda 3

This page was generated from /home/docs/checkouts/readthedocs.org/user_builds/blm/checkouts/latest/docs/tasks/task1.ipynb. Interactive online version: Slideshow:

Understanding Surface Parameters (1): Albedo

Objectives

1. Review Python programming.

2. Prepare plots of surface radiation and energy balance using assigned AmeriFlux observations under two scenarios: one clear day and one cloudy day.

Tips:

1. Check the fluxes make sense.

2. Think about the time of day.

3. Calculate the albedo through the day for both days.

4. How does the albedo vary through the day?

5. How do your site values compare to the literature?

6. How does the albedo vary through the year?

7. If you had to use one albedo value to model the site, what would it be?

8. For scientific background see Albedo and Reading List

Tasks

load necessary packages

# These packages are necessary later on. Load all the packages in one place for consistency
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from pathlib import Path

load data

#The path of the directory where all AMF data are
#path_dir = Path.cwd()/'data'/'AMF-clean'
path_dir = Path.cwd()/'data'/'1'

name_of_site = 'CA-Obs_clean.csv.gz'
path_data = path_dir/name_of_site
path_data.resolve()

PosixPath('/Users/sunt05/Dropbox/6.Repos/BLM/data/1/CA-Obs_clean.csv.gz')

df_data = pd.read_csv(path_data, index_col='time',parse_dates=['time'])

---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
<ipython-input-4-010038206d98> in <module>
----> 1 df_data = pd.read_csv(path_data, index_col='time',parse_dates=['time'])

/anaconda3/lib/python3.7/site-packages/pandas/io/parsers.py in parser_f(filepath_or_buffer, sep, delimiter, header, names, index_col, usecols, squeeze, prefix, mangle_dupe_cols, dtype, engine, converters, true_values, false_values, skipinitialspace, skiprows, skipfooter, nrows, na_values, keep_default_na, na_filter, verbose, skip_blank_lines, parse_dates, infer_datetime_format, keep_date_col, date_parser, dayfirst, cache_dates, iterator, chunksize, compression, thousands, decimal, lineterminator, quotechar, quoting, doublequote, escapechar, comment, encoding, dialect, error_bad_lines, warn_bad_lines, delim_whitespace, low_memory, memory_map, float_precision)
    683         )
    684
--> 685         return _read(filepath_or_buffer, kwds)
    686
    687     parser_f.__name__ = name

/anaconda3/lib/python3.7/site-packages/pandas/io/parsers.py in _read(filepath_or_buffer, kwds)
    455
    456     # Create the parser.
--> 457     parser = TextFileReader(fp_or_buf, **kwds)
    458
    459     if chunksize or iterator:

/anaconda3/lib/python3.7/site-packages/pandas/io/parsers.py in __init__(self, f, engine, **kwds)
    893             self.options["has_index_names"] = kwds["has_index_names"]
    894
--> 895         self._make_engine(self.engine)
    896
    897     def close(self):

/anaconda3/lib/python3.7/site-packages/pandas/io/parsers.py in _make_engine(self, engine)
   1133     def _make_engine(self, engine="c"):
   1134         if engine == "c":
-> 1135             self._engine = CParserWrapper(self.f, **self.options)
   1136         else:
   1137             if engine == "python":

/anaconda3/lib/python3.7/site-packages/pandas/io/parsers.py in __init__(self, src, **kwds)
   1915         kwds["usecols"] = self.usecols
   1916
-> 1917         self._reader = parsers.TextReader(src, **kwds)
   1918         self.unnamed_cols = self._reader.unnamed_cols
   1919

pandas/_libs/parsers.pyx in pandas._libs.parsers.TextReader.__cinit__()

pandas/_libs/parsers.pyx in pandas._libs.parsers.TextReader._setup_parser_source()

/anaconda3/lib/python3.7/gzip.py in __init__(self, filename, mode, compresslevel, fileobj, mtime)
    161             mode += 'b'
    162         if fileobj is None:
--> 163             fileobj = self.myfileobj = builtins.open(filename, mode or 'rb')
    164         if filename is None:
    165             filename = getattr(fileobj, 'name', '')

FileNotFoundError: [Errno 2] No such file or directory: '/Users/sunt05/Dropbox/6.Repos/BLM/data/1/CA-Obs_clean.csv.gz'

Let’s look at the site info

all_sites_info =  pd.read_csv('data/site_info.csv')
site_info=all_sites_info[all_sites_info['Site Id'] == name_of_site.split('_clean')[0]]
site_info

examine data

Loading the first 5 rows of the dataframe

df_data.head(5)

Let’s print some statistics of the data

For names used by AmeriFlux for variables see AmeriFlux Key variables for analysis

df_data.describe()

df_data.loc[:,['SWIN','LWIN','SWOUT','LWOUT']].plot(figsize=(15,4))

df_data.loc[:,['LE','H']].plot(figsize=(15,4))

df_data.loc[:,['TA']].plot(figsize=(15,4))
df_data.loc[:,['P']].plot(figsize=(15,4))
df_data.loc[:,['WS']].plot(figsize=(15,4))
df_data.loc[:,['RH']].plot(figsize=(15,4))
df_data.loc[:,['PA']].plot(figsize=(15,4))

plot data

radiation balance

date = '2001 10 21'
df_data.loc[:,['SWIN','LWIN','SWOUT','LWOUT']].dropna().loc[date].plot()

examining if NETRAD == (SWIN-SWOUT)+(LWIN-LWOUT)

NetSW = df_data.loc[date,'SWIN']-df_data.loc[date,'SWOUT']
NetLW = df_data.loc[date,'LWIN']-df_data.loc[date,'LWOUT']
Netrad_calc=NetSW+NetLW
df_data.loc[date,'NETRAD'].plot(color='r',label='NETRAD from data')
Netrad_calc.plot(color='b',marker='o',label='NETRAD calculated')
plt.legend()

surface energy balance

df_data.loc[:,['NETRAD','H','LE']].dropna().loc[date].plot()

calculate albedo

surface albedo

ser_alb=df_data['SWOUT']/df_data['SWIN']

ser_alb[ser_alb.between(0,1)&(df_data['SWIN']>5)].loc[date].plot(marker='o')

Note

Please report issues or ask questions about this site on the GitHub page.

Visit to Observatory

Objective

The purpose of this visit is to think about how instrumentation is deployed and its fetch. See more details here: URAO

Tip

1. Read Safety information about visiting the URAO

2. Read the Risk Assessment for visiting the Observatory (on BB).

3. Sign the Risk Assessment for visiting the Observatory.

Activity in URAO

1. Identify the instrumentation to directly measure all terms of the surface energy balance, radiation balance, temperatures, wind speed, wind direction, boundary layer height and/or cloud height and various land surface model parameters. What are the instrument details, including name, manufacturer, sample rate, averaging time, height, orientation (e.g. wind direction of sonic anemometer), other variables observed or derived by/from an instrument

2. Look at the surroundings of the University of Reading Atmospheric Observatory (URAO). Identify in what wind directions will the sensors have different fetch (nature of upwind surface).

Activity after URAO visit

1. For your site in the AmeriFlux network, access key papers about your site to determine how the fluxes and other variables are observed. Create a summary table that includes:

   • Instrument details, including name, manufacturer, sample rate, averaging time, height, orientation (e.g. wind direction of sonic anemometer)

   • Photograph (from web) of the instruments (i.e. model,etc) you can work in teams to find these.

   • General method used to obtain the flux

2. Using Google/Bing Maps (or other equivalent) to identify what upwind conditions are relative to the orientation of the instruments. Think carefully about which sensors this is critical for.

   • Create an annotated figure for your site with sectors identified.

3. Prepare a document that summarises the information above for the start of next class.

   • The document should have figure and table captions (if you are unclear review the journal papers you are reading about your site).

   • Note: A figure/table caption should be sufficient for a reader to interpret the material without reading the body of the text.

This page was generated from /home/docs/checkouts/readthedocs.org/user_builds/blm/checkouts/latest/docs/tasks/task2-1.ipynb. Interactive online version: Slideshow:

Understanding Surface Properties (2): LAI

Objectives

1. Understand the dynamics of LAI at different temporal scales (e.g., sub-daily, seasonal)

2. Understand the effects of LAI dynamics on albedo

Tasks

load necessary packages

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from datetime import datetime
from pathlib import Path
import datetime

Name of the site and group number

loading data for albedo calculation

group_number=2
path_dir = Path.cwd()/'data'/f'{group_number}'
# examine available files in your folder
list(path_dir.glob('*gz'))

[PosixPath('/Users/hamidrezaomidvar/Desktop/BLM/docs/tasks/data/2/US-MMS_clean.csv.gz')]

name_of_site='US-MMS' #

# load dataset
site_file=name_of_site+'_clean.csv.gz'
path_data = path_dir/site_file
df_data = pd.read_csv(path_data, index_col='time',parse_dates=['time'])
df_data.head()

	WS	RH	TA	PA	WD	P	SWIN	LWIN	SWOUT	LWOUT	NETRAD	H	LE
time
1999-01-01 01:00:00	2.98	NaN	NaN	99.0	95.50	0.0	0.0	190.0	0.0	250.7	-60.7	-21.310	-0.350
1999-01-01 02:00:00	2.93	NaN	NaN	99.0	108.71	0.0	0.0	188.0	0.0	248.8	-60.8	-14.473	-0.321
1999-01-01 03:00:00	2.96	NaN	NaN	99.0	122.46	0.0	0.0	187.5	0.0	247.6	-60.1	-16.784	-0.786
1999-01-01 04:00:00	3.03	NaN	NaN	99.1	118.99	0.0	0.0	187.6	0.0	247.7	-60.1	-12.367	-0.191
1999-01-01 05:00:00	3.63	NaN	NaN	99.2	124.80	0.0	0.0	186.4	0.0	246.9	-60.5	-23.069	-0.650

Take a look at the raditation data

fontsize=15
df_data.loc[:,['SWIN','LWIN','SWOUT','LWOUT']].plot(figsize=(12,4),fontsize=fontsize)
plt.ylabel('Flux (W m$^{-2}$)',fontsize=fontsize)
plt.xlabel('Time',fontsize=fontsize)

Text(0.5, 0, 'Time')

Loading LAI data + some data cleaning

def DOY_to_datetime(row):
    year=int(row['modis_date'][1:5])
    DOY=int(row['modis_date'][5:])
    return datetime.datetime(year, 1, 1) + datetime.timedelta(DOY - 1)

df_LAI=pd.read_csv('data/MODIS_LAI_AmeriFlux/statistics_Lai_500m-'+name_of_site+'.csv')
df_LAI.columns=['product']+[i.split(' ')[1] for i in df_LAI.columns if i!='product']
df_LAI=df_LAI.filter(['modis_date','value_mean'])


df_LAI.set_index(df_LAI.apply(DOY_to_datetime,axis=1),inplace=True)
df_LAI.drop('modis_date',axis=1,inplace=True)
df_LAI.head()

	value_mean
2002-07-04	5.1042
2002-07-08	3.2211
2002-07-12	3.6012
2002-07-16	4.2047
2002-07-20	4.0568

fontsize=15
df_LAI.plot(legend=False,figsize=(12,4),fontsize=fontsize)
plt.ylabel('LAI',fontsize=fontsize)
plt.xlabel('Time',fontsize=fontsize)

Text(0.5, 0, 'Time')

examine results

Choose the desire period

start_period='2006-01-01'
end_period='2007-01-01'

df_data_period=df_data.loc[start_period:end_period]
df_LAI_period=df_LAI.loc[start_period:end_period]

Alebdo calculation

ser_alb=df_data_period['SWOUT']/df_data_period['SWIN']
ser_alb=ser_alb[ser_alb.index.time==datetime.time(12, 0)]
ser_alb_clean=ser_alb[ser_alb.between(0,1)&(df_data_period['SWIN']>5)]


plt.rcParams.update({'font.size': 15})
plt.figure(figsize=(12,4))
plt.scatter(ser_alb_clean.index,ser_alb_clean)
plt.ylabel('Albedo')
plt.xlabel('Time')

Text(0.5, 0, 'Time')

explore the relationship between albedo and LAI

plt.rcParams.update({'font.size': 15})
fig,axs=plt.subplots(2,1,figsize=(12,8))
fig.subplots_adjust(hspace=0)
ax1=axs[0]
ax2=axs[1]

ax1.scatter(ser_alb_clean.index,ser_alb_clean)
ax1.set_ylabel('Albedo')
ax1.set_xlim([start_period,end_period])
ax1.set_xticks([])
ax1.set_title(name_of_site)


ax2.plot(df_LAI_period.index,df_LAI_period)
ax2.set_ylabel('LAI')
ax2.set_xlabel('Time')
ax2.set_xlim([start_period,end_period])

(732312.0, 732677.0)

This page was generated from /home/docs/checkouts/readthedocs.org/user_builds/blm/checkouts/latest/docs/tasks/task3.ipynb. Interactive online version: Slideshow:

Understanding Surface Properties (3): Surface Roughness

Objectives

1. Understand how to quantify surface roughness

2. Understand variability of surface roughness

Tips:

1. Eddy Covariance observations allow several important parameters to be determined:

• Zero plane displacement length \(d\) can be calculated based on vegetation height (think about changes in LAI and porosity). If there is a profile of anemometers (e.g. as at the Observatory) then this can be calculated using that.

• Roughness length \(z_0\) can be calculated from the EC data (EddyCovariance) by rearranging the logarithmic law from the neutral wind profile WindProfile).

2. Do these vary with season? Wind Direction? Create a graph to analyse the data (think about wind directions).

3. How do these values compare to the literature?

4. Calculate the aerodynamic resistance \(r_a\).

5. For scientific background see Roughness and Reading List

Tasks

load necessary packages

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from pathlib import Path

loading data

group_number = 10
path_dir = Path.cwd() / 'data' / f'{group_number}'
# examine available files in your folder
list(path_dir.glob('*gz'))

[PosixPath('/Users/sunt05/Dropbox/6.Repos/BLM-task3/data/10/US-Oho_clean.csv.gz'),
 PosixPath('/Users/sunt05/Dropbox/6.Repos/BLM-task3/data/10/CA-TPD_clean.csv.gz')]

# specify the site name
name_of_site = 'CA-TPD'

# load measurement heights
df_zmeas = pd.read_csv('data/measurement_height.csv', index_col=0)
df_zmeas.loc[name_of_site].iloc[:]

	Variable	Start_Date	Height	Instrument_Model	Instrument_Model2	Comment	BASE_Version
Site_ID
CA-TPD	CO2_1_1_1	2012.0	35.70	GA_CP-LI-COR LI-7200	NaN	LI-7200 IRGA)	02-May
CA-TPD	CO2_1_2_1	2012.0	16.00	GA_SR-LI-COR LI-820	NaN	LI820_CO2_Cnpy_16m (ppm); CO2	02-May
CA-TPD	FC	2012.0	35.70	GA_CP-LI-COR LI-7200	NaN	NaN	02-May
CA-TPD	FH2O	2012.0	NaN	NaN	NaN	NaN	02-May
CA-TPD	LE	2012.0	35.70	GA_CP-LI-COR LI-7200	NaN	NaN	02-May
CA-TPD	PPFD_IN_PI_F_1_1_1	2012.0	35.13	RAD-PAR Quantum	NaN	NaN	02-May
CA-TPD	P_CUM	2012.0	2.00	NaN	NaN	NaN	02-May
CA-TPD	PA_PI_F	2012.0	2.00	NaN	NaN	NaN	02-May
CA-TPD	G_1_1_1	2012.0	0.00	SOIL_H-Plate	NaN	NaN	02-May
CA-TPD	G_2_2_1	2012.0	-0.03	SOIL_H-Plate	NaN	NaN	02-May
CA-TPD	G_3_2_1	2012.0	-0.03	SOIL_H-Plate	NaN	NaN	02-May
CA-TPD	G_4_2_1	2012.0	-0.03	SOIL_H-Plate	NaN	NaN	02-May
CA-TPD	G_1_2_1	2012.0	-0.03	SOIL_H-Plate	NaN	NaN	02-May
CA-TPD	GPP_PI_F	2012.0	35.70	GA_CP-LI-COR LI-7200	NaN	NaN	02-May
CA-TPD	H	2012.0	35.70	GA_CP-LI-COR LI-7200	NaN	NaN	02-May
CA-TPD	H_PI_F	2012.0	35.70	GA_CP-LI-COR LI-7200	NaN	NaN	02-May
CA-TPD	LE_PI_F	2012.0	35.70	GA_CP-LI-COR LI-7200	NaN	NaN	02-May
CA-TPD	NEE_PI	2012.0	35.70	GA_CP-LI-COR LI-7200	NaN	NaN	02-May
CA-TPD	NEE_PI_F	2012.0	35.70	GA_CP-LI-COR LI-7200	NaN	NaN	02-May
CA-TPD	PPFD_IN_1_1_1	2012.0	35.13	RAD-PAR Quantum	NaN	Kipp and Zonen	02-May
CA-TPD	PPFD_IN_1_2_1	2012.0	2.00	RAD-PAR Quantum	NaN	NaN	02-May
CA-TPD	PPFD_OUT	2012.0	35.13	RAD-PAR Quantum	NaN	Kipp and Zonen	02-May
CA-TPD	P	2012.0	2.00	PREC-TipBucGauge	NaN	CSI (heated and Alter Shielded) in a forest op...	02-May
CA-TPD	PA	2012.0	2.00	PRES-ElectBar	NaN	Young Co	02-May
CA-TPD	RECO_PI_F	2012.0	35.70	GA_CP-LI-COR LI-7200	NaN	NaN	02-May
CA-TPD	SW_IN	2012.0	35.50	RAD-Pyrrad-SW+LW	NaN	Kipp & Zonen CNR4	02-May
CA-TPD	SW_IN_PI_F	2012.0	35.50	RAD-Pyrrad-SW+LW	NaN	CNR1_GlobalShortwaveRad_AbvCnpy_36m (W/m2); Ga...	02-May
CA-TPD	LW_IN	2012.0	35.50	RAD-Pyrrad-SW+LW	NaN	Kipp & Zonen CNR4 Net Radiometer	02-May
CA-TPD	LW_OUT	2012.0	35.50	RAD-Pyrrad-SW+LW	NaN	Kipp & Zonen CNR4 Net Radiometer	02-May
CA-TPD	SW_OUT	2012.0	35.50	RAD-Pyrrad-SW+LW	NaN	Kipp & Zonen CNR4 Net Radiometer	02-May
...	...	...	...	...	...	...	...
CA-TPD	SWC_2_1_1	2012.0	-0.02	SWC-TDR	NaN	NaN	02-May
CA-TPD	SWC_2_3_1	2012.0	-0.10	SWC-TDR	NaN	NaN	02-May
CA-TPD	SWC_2_6_1	2012.0	-1.00	SWC-TDR	NaN	NaN	02-May
CA-TPD	SWC_2_4_1	2012.0	-0.20	SWC-TDR	NaN	NaN	02-May
CA-TPD	SWC_2_2_1	2012.0	-0.05	SWC-TDR	NaN	NaN	02-May
CA-TPD	SWC_2_5_1	2012.0	-0.50	SWC-TDR	NaN	NaN	02-May
CA-TPD	SWC_1_3_1	2012.0	-0.10	SWC-TDR	NaN	NaN	02-May
CA-TPD	SWC_1_6_1	2012.0	-1.00	SWC-TDR	NaN	NaN	02-May
CA-TPD	SWC_1_4_1	2012.0	-0.20	SWC-TDR	NaN	NaN	02-May
CA-TPD	SWC_1_2_1	2012.0	-0.05	SWC-TDR	NaN	NaN	02-May
CA-TPD	SWC_1_5_1	2012.0	-0.50	SWC-TDR	NaN	NaN	02-May
CA-TPD	TS_2_3_1	2012.0	-0.10	TEMP-ElectResis	NaN	NaN	02-May
CA-TPD	TS_2_6_1	2012.0	-1.00	TEMP-ElectResis	NaN	NaN	02-May
CA-TPD	TS_2_4_1	2012.0	-0.20	TEMP-ElectResis	NaN	NaN	02-May
CA-TPD	TA	2012.0	36.60	TEMP-ElectResis	NaN	Campbell Scientific Inc (CSI)	02-May
CA-TPD	TS_2_2_1	2012.0	-0.05	TEMP-ElectResis	NaN	NaN	02-May
CA-TPD	TS_2_5_1	2012.0	-0.50	TEMP-ElectResis	NaN	NaN	02-May
CA-TPD	TS_1_3_1	2012.0	-0.10	TEMP-ElectResis	NaN	NaN	02-May
CA-TPD	TS_1_6_1	2012.0	-1.00	TEMP-ElectResis	NaN	NaN	02-May
CA-TPD	TS_1_4_1	2012.0	-0.20	TEMP-ElectResis	NaN	NaN	02-May
CA-TPD	TS_1_2_1	2012.0	-0.05	TEMP-ElectResis	NaN	Campbell Scientific Inc (CSI)	02-May
CA-TPD	TS_1_5_1	2012.0	-0.50	TEMP-ElectResis	NaN	NaN	02-May
CA-TPD	TS_PI_F_1_1_1	2012.0	-0.02	TEMP-ElectResis	NaN	NaN	02-May
CA-TPD	TS_PI_F_1_2_1	2012.0	-0.05	TEMP-ElectResis	NaN	NaN	02-May
CA-TPD	USTAR	2012.0	35.70	SA-Campbell CSAT-3	NaN	NaN	02-May
CA-TPD	VPD_PI	2012.0	35.30	RH-Capac	NaN	NaN	02-May
CA-TPD	VPD_PI_F	2012.0	35.30	RH-Capac	NaN	NaN	02-May
CA-TPD	WD	2012.0	36.60	WIND-VaneAn	NaN	NaN	02-May
CA-TPD	WS	2012.0	36.60	WIND-VaneAn	NaN	NaN	02-May
CA-TPD	WS_PI_F	2012.0	36.60	WIND-VaneAn	NaN	NaN	02-May

73 rows × 7 columns

from the table above, we know the measurement height for wind speed WS is 36.6 m.

We will use the value in the following calculation.

# write down measurement height of wind speed
z_meas = 36.6

# load dataset
site_file = name_of_site + '_clean.csv.gz'
path_data = path_dir / site_file
df_data = pd.read_csv(path_data, index_col='time', parse_dates=['time'])
df_data.head()

	WS	RH	TA	PA	WD	P	SWIN	LWIN	SWOUT	LWOUT	NETRAD	H	LE	USTAR	ZL
time
2012-01-01 11:00:00	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	-86.133	12.683	1.170	NaN
2012-01-01 11:30:00	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	-71.325	28.933	1.188	NaN
2012-01-01 12:00:00	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	-64.607	29.349	1.188	NaN
2012-01-01 12:30:00	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	-94.362	36.909	1.295	NaN
2012-01-01 13:00:00	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	-74.555	51.049	1.212	NaN

stability correction

You need to calculate Obukhov Stavility parameter as following:

\[L=-\frac{u_{*}^{3}}{k \frac{g}{T} \frac{Q_{H}}{\rho c_{p}}}\]

\[\text{Obukhov stability parameter }=\frac{z-d}{L}\]

where \(u_*\) is friction velocity; \(g\) is gravity; \(T\) is air temperature; \(k\) is von Kármán’s ‘constant’ (0.4); \(Q_H\) is turbulent sensible heat flux; \(\rho\) is air density; \(z\) is measurement height; and \(d\) is zero desplacement height

We use different functions for this. The functions are already implemented. You can use these functions in the future for different part of the model. Here some example how to use them.

from utility import cal_vap_sat, cal_dens_dry, cal_dens_vap, cal_cpa, cal_dens_air, cal_Lob

# saturation vapour pressure [hPa]
cal_vap_sat(0.0003, 1013)

6.118926455724389

# density of dry air [kg m-3]
cal_dens_dry(98, 23, 1013)

1.1591455092830807

# density of vapour [kg m-3]
cal_dens_vap(98, 23, 1013)

0.020167600169459326

# specific heat capacity of air mass [J kg-1 K-1]
cal_cpa(23, 30, 1013)

1010.1394406405146

# air density [kg m-3]
cal_dens_air(1013, 32)

1.1564283022625206

# Obukhov length
cal_Lob(500, 0.0001, 23, 90, 1020)

-0.00018481221694860222

Now we can calculate Obukhov length (\(L\)). Let’s integrate all together:

Here you need to specify the height of the vegetation h_sfc . You can find this from the paper related to your site. Note that the height should be smaller than the measurement height

# select valid values
df_val = df_data.loc[:, ['H', 'USTAR', 'TA', 'RH', 'PA', 'WS']].dropna()

# calculate Obukhov length
ser_Lob = df_val.apply(
    lambda ser: cal_Lob(ser.H, ser.USTAR, ser.TA, ser.RH, ser.PA * 10), axis=1)

# zero-plane displacement: estimated using rule f thumb `d=0.7*h_sfc`
h_sfc = 10
z_d = 0.7 * h_sfc

if z_d >= z_meas:
    print(
        'vegetation height is greater than measuring height. Please fix this before continuing'
    )

# calculate stability scale
ser_zL = (z_meas - z_d) / ser_Lob

# determine periods under quasi-neutral conditions
limit_neutral = 0.01
ind_neutral = ser_zL.between(-limit_neutral, limit_neutral)
df_val.loc[ind_neutral, 'USTAR'].plot.hist()

<matplotlib.axes._subplots.AxesSubplot at 0x1158fd940>

calculate roughness length \(z_0\) and zero plane displacement \(d\)

theoretical basis

We will employ the followng equation in our calculations:

\[\bar{u}(z)=\frac{u_{*}}{k} \ln \left[\frac{(z-d)}{z_{0}}\right]\]

where \(\bar{u}(z)\) is the measured wind speed at a height \(z\), \(u_{*}\) the friction velocity, \(d\) the zero-plane displacement, \(k\) the von Karman’s constant (a value of 0.4 used in the following calculation).

calculation

We will use scipy.optimize.curve_fit to determine the parameters \(z_0\) and \(d\) in the model above.

from scipy.optimize import curve_fit

To use curve_fit, we need to define a function repensenting the model above:

def func_uz(ustar, z0):
    z = z_meas
    d = 0.7 * h_sfc
    k = 0.4
    uz = ustar / k * np.log((z - d) / z0)
    return uz

# choose the necessary variables
df_sel = df_val.loc[ind_neutral, ['WS', 'USTAR']].dropna()

# set up input data for curve fitting
ser_ustar = df_sel.USTAR
ser_ws = df_sel.WS
# use `curve_fit` to get parameters of interest
z0_fit, _ = curve_fit(func_uz, ser_ustar, ser_ws, bounds=[(0), (np.inf)])

# examine the results
res_fit = func_uz(ser_ustar, z0_fit).rename('sim')
res_obs = ser_ws.copy().rename('obs')
df_comp = pd.concat([res_obs, res_fit], axis=1)
ax = df_comp.plot.scatter(x='obs', y='sim', label=f'$z_0=${z0_fit[0]:.2f}')
lim = ax.get_xlim()
ax.plot(lim, lim, color='r', linestyle='--')
ax.set_aspect('equal')

[<matplotlib.lines.Line2D at 0x127a4a048>]

Calculate aerodynamic resistance \(r_a\)

write down the equation used for this calculation

equation with explanations of symbol:

\[equation\ here\]

calculation

# calculation code

examination of variability

# calculation code

This page was generated from /home/docs/checkouts/readthedocs.org/user_builds/blm/checkouts/latest/docs/tasks/task4.ipynb. Interactive online version: Slideshow:

Understanding Surface Properties (4): Heat Storage and Surface Resistance

Objectives

1. Understand how to calculate ground heat flux.

2. Understand how to back-calculate surface resistance.

Calculation of heat storage: OHM (Objective Hysteresis Model) and a simple scheme

The objective hysteresis model helps us estimate the ground surface heat flux using the measured net radiation data. The equation for this model is as follows:

\[\Delta Q_{S}=a_{1} Q^{*}+a_{2} \frac{\partial Q^{*}}{\partial t}+a_{3}\]

where \(\Delta Q_S\) is the ground heat flux (W m\(^{-2}\)), \(Q^*\) is the net radiation (W m\(^{-2}\)), and \(a_{1-3}\) are the coefficients to be determined from measurements

In this task, we focus on a simpler model instead of the above question using just the first term:

\[\Delta Q_{S}=a_{1} Q^{*}\]

The aim is to find \(a_1\) using the measurements.

Tasks

Load necessary packages

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from pathlib import Path

Loading data

group_number = 7
path_dir = Path.cwd() / 'data' / f'{group_number}'
# examine available files in your folder
list(path_dir.glob('*gz'))

[PosixPath('/Users/hamidrezaomidvar/Desktop/BLM-task4/data/7/US-Whs_clean.csv.gz'),
 PosixPath('/Users/hamidrezaomidvar/Desktop/BLM-task4/data/7/US-NC1_clean.csv.gz')]

# specify the site name
name_of_site = 'US-Whs'

# load dataset
site_file = name_of_site + '_clean.csv.gz'
path_data = path_dir / site_file
df_data = pd.read_csv(path_data, index_col='time', parse_dates=['time'])
df_data.head()

	WS	RH	TA	PA	WD	P	SWIN	LWIN	SWOUT	LWOUT	NETRAD	H	LE	USTAR	ZL
time
2007-06-29 13:30:00	4.063	12.6475	34.790	86.3667	273.75	0.0	1096.50	NaN	NaN	NaN	609.9	339.3645	12.2582	0.4036	NaN
2007-06-29 14:00:00	3.636	12.2413	35.140	86.3323	322.30	0.0	1057.00	NaN	NaN	NaN	562.8	290.9766	14.5989	0.5052	NaN
2007-06-29 14:30:00	3.107	12.0811	35.490	86.3048	316.75	0.0	1107.50	NaN	NaN	NaN	596.4	359.2087	16.1474	0.3566	NaN
2007-06-29 15:00:00	4.302	11.6801	36.190	86.2589	292.30	0.0	1108.00	NaN	NaN	NaN	613.5	360.5630	10.0565	0.5340	NaN
2007-06-29 15:30:00	3.676	13.7432	34.785	86.2429	209.20	0.0	496.15	NaN	NaN	NaN	232.0	200.0689	64.8520	0.4503	NaN

Assuming surface energy closure in measurements is perfect, i.e., \(Q^*=Q_H+Q_E+\Delta Q_S\)

Note: conduct the derivation at different sites separately so you can get two sets of coefficients for later comparison.

Calculating \(\Delta Q_S\)

df=df_data.filter(['NETRAD','H','LE'])
df.dropna(inplace=True)
QSTAR=df.NETRAD
QH=df.H
QE=df.LE

QS=QSTAR-QH-QE


QS.plot(figsize=(15,4))
plt.ylabel('$\Delta Q_S$')

Text(0, 0.5, '$\\Delta Q_S$')

Finding \(a_1\)

from scipy.stats import linregress

fig, ax = plt.subplots(1, 1)
slope = linregress(QSTAR, QS).slope
intercept=linregress(QSTAR, QS).intercept

ax.scatter(QSTAR, QS)
lim = ax.get_xlim()
plt.plot(lim, [x * slope+intercept for x in lim], color='r', linestyle='--')
plt.title(f'$a_1 = {slope:.2f}$')
plt.ylabel('$\Delta Q_S$')
plt.xlabel('$Q^*$')

Text(0.5, 0, '$Q^*$')

Compare derived OHM coefficients between sites

# your code here

Calculate \(\Delta Q_S\) using provided OHM coefficients

Use the values provided in this page of SUEWS manual

# your code here

examine the SEB closure using calculated

Think about proper metrics for quantitative examination.

# your code here

calculate surface resistances \(r_s\)

\(r_s\) can be calculated as following

\[r_{s}=\left(\frac{s \beta}{\gamma}-1\right) r_{a}+\frac{C_{a} V}{\gamma Q_{E}}\]

where

\[\beta=\frac{Q_{H}}{Q_{E}}\]

\[C_{a}=\rho c_{p}\]

\[V=e_{s}-e_{a}\]

\[\gamma = \text{psychrometric constant}\]

\[s=\text{slope of saturation evaporation pressure curve vs air temperature}\]

We are going to calculate each part separately and eventually integrate all parts together to get \(r_s\)

Calculate \(r_a\)

You should calculate \(r_a\) using task 3 and the function that you implemented before for the site: similarly for \(z_0\) and \(d\).

# your code here

Calculating \(\beta=\frac{Q_{H}}{Q_{E}}\)

df_val = df_data.loc[:, ['LE','H', 'USTAR', 'TA', 'RH', 'PA', 'WS']].dropna()

# your code here

Calculating \(C_{a}=\rho c_{p}\)

For this, you can use the available functions in utility.py

Important:

In these functions, pressure is in hPa; so you need to convert the pressure from kPa to hPa

df_val.PA *=10

from utility import cal_dens_air, cal_cpa

#rho = cal_dens_air(Press_hPa, Temp_C)  #use this to calculate rho
#cp = cal_cpa(Temp_C, RH_pct, Press_hPa) #use this to calculate cp

Calculating \(s\)

# This package is needed
!pip install atmosp

from utility import cal_des_dta
# s = cal_des_dta(Temp_C, Press_hpa) #use this to calculate s

Calculating \(V=e_{s}-e_{a}\)

from utility import cal_vpd
# V = cal_vpd(Temp_C, RH_pct, Press_hPa)  #use this to calculate V

Calculating \(\gamma\)

from utility import cal_gamma
#gamma = cal_gamma(Press_hPa) #use this to calculate gamma

Integrate: calculate \(r_s\)

# your code here

This page was generated from /home/docs/checkouts/readthedocs.org/user_builds/blm/checkouts/latest/docs/tasks/task5.ipynb. Interactive online version: Slideshow:

Modelling \(Q_E\) and \(Q_H\) with Penman-Monteith

Objectives

1. Calculate \(Q_E\) and \(Q_H\) using Penman-Monteith

2. Evaluate the above results with EC measurements

Tips:

1. The Penman Monteith equation assumes energy balance closure at each model time step.

2. Identify parameter values used initially.

3. Use the calculated roughness length and the canopy resistances to calculate the PM fluxes (QH and QE) for the same period. What is the impact on the evaluation statistics?

Task explained

basic equations

The Penman-Monteith equation:

\[Q_{E} = \frac{s\left( Q^{*} - Q_{G} \right) + \rho c_{p}(e_{s} - e)/r_{H}} {s + \gamma\ (1 + r_{s}/r_{a})}\]

where

\(s\): the slope of saturation vapour pressure curve at temperature \(T_a\);

\(Q^*\): net all-wave radiation;

\(Q_G\): ground heat flux;

\(\rho\): air density;

\(c_p\): heat capacity of air;

\(e_s\): saturation vapour pressure at temperature \(T_a\);

\(e\): actual vapour pressure;

\(\gamma\): psychrometric constant;

\(r_s\): surface resistance;

\(r_H\): aerodynamic resistance for heat transfer;

\(r_a\): aerodynamic resistance for momentum transfer.

The surface energy balance equation:

\[Q^*-Q_G=Q_E+Q_H\]

where \(Q^*\) is the net all-wave radiation, \(Q_G\) the ground heat flux, \(Q_E\) the latent heat flux and \(Q_H\) the sensible heat flux.

specific aims

Before starting coding, we need to interpret the problem into the following major task: With atmospheric forcing conditions and other variables in your specific AmeriFlux dataset, one needs

1. to model \(Q_E\) using the above equation and \(Q_H\) using surface energy balance (SEB) equation； and

2. to evaluate the modelled \(Q_E\) and \(Q_H\) against their observed counterparts using quantitative metrics.

assumptions

For simplicity, we make the following assumptions:

1. perfect SEB closure;

2. known outgoing longwave radiation; and

3. equivalence between \(r_H\) and \(r_a\): i.e. :math:`r_H=r_a`.

guide questions for modelling \(Q_E\)

As \(Q_H\) would be easily obtained as residual of other terms (i.e., \(Q_H=Q^*-Q_G-Q_E\)), let’s first focus on how to get \(Q_E\) by thinking the following questions:

1. What are the known and unknown terms in the above given the progress so far?

2. For each of the unknown terms, is it related to atmospheric forcing or surface properties?

input data we have

As a reminder, we know that we have the following atmospheric forcing variables available in the AmeriFlux dataset:

1. incoming solar radiation \(K_\downarrow\);

2. incoming longwave radiation \(L_\downarrow\);

3. outgoing longwave radiation \(L_\uparrow\);

4. air temperature \(T_a\);

5. relative humidity \(RH\);

6. atmospheric pressure \(p\);

7. wind speed \(u\); and

8. precipitation \(P\) (note this is NOT used in this task).

Other variables are provided either as auxiliary data to derive surface properties (e.g., friction velocity \(u_*\)) or as observed values for evaluation of modeled variables (e.g., sensible heat flux \(Q_H\)).

Now let’s move onto the calculations below: more specific hints are provided in comments of cells below.

Tasks

load necessary packages

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from pathlib import Path

Note:

The following cell enables the autoreload module, which allows auto-reloading of external packages. This would be particularly useful if we make changes in the utility.py file: such changes will take immediate effec in this notebook without restarting the Jupyter kernal as one would usually do.

%load_ext autoreload
%autoreload 2

load data

group_number = 5
path_dir = Path.cwd() / 'data' / f'{group_number}'
# examine available files in your folder
list(path_dir.glob('*gz'))

[PosixPath('/Users/sunt05/Dropbox/6.Repos/task-5-blm2019/data/5/US-Slt_clean.csv.gz'),
 PosixPath('/Users/sunt05/Dropbox/6.Repos/task-5-blm2019/data/5/US-UMB_clean.csv.gz')]

# specify the site name
name_of_site = 'US-UMB'

# load dataset
site_file = name_of_site + '_clean.csv.gz'
path_data = path_dir / site_file
df_data = pd.read_csv(path_data, index_col='time', parse_dates=['time'])
df_data.head()

Calculate \(Q_E\) using Penman-Monteith

calculate \(Q^*\) using surface radiation balance

Note: we’ve got albedo in Task 1.

# your code

calculate \(Q_G\) using a simple linear model

Note: we’ve done this in Task 4.

# your code

calculate \(r_a\)

Note: we’ve done this in Task 3.

# your code

calculate \(r_s\)

Note: we’ve done this in Task 4.

# your code

calculate other atmospheric states

Note: related functions are all provided in the utility.py file.

# slope `s`:
# s = cal_des_dta(Temp_C, Press_hPa)

# psychrometric constant `gamma`:
# gamma = cal_gamma(Press_hPa)

# air density `rho`
# rho = cal_dens_air(Press_hPa, Temp_C)

# specific heat capacity of air mass `cp`
# cp = cal_cpa(Temp_C, RH_pct, Press_hPa)

# vapour pressure deficit `es-e`, or `vpd`
# vpd = cal_vpd(Temp_C, RH_pct, Press_hPa)

integrate the above results to get \(Q_E\)

# your code

calculate \(Q_H\) using SEB equation

# your code

evaluation of modelled \(Q_E\) and \(Q_H\)

# calculate MAE to see the overal performance

# calculate MBE to see if your model is overestimating or underestimating the variables

# plot observed and modelled values on the same chart with a 1:1 line
# note: we've done this in task 3 for `ra`

This page was generated from /home/docs/checkouts/readthedocs.org/user_builds/blm/checkouts/latest/docs/tasks/task6_extra-MT49E.ipynb. Interactive online version: Slideshow:

MT49E Extra: Sensitivity Analysis of Modelled \(Q_E\)

Objectives

1. To assess which parameters and varaiables influence your model most

Tips:

1. What are the effects of systematically perturbing \(Q^*\), \(Q_G\), canopy resistance \(r_c\), roughness length \(z_0\) in the Penman Monteith equation?

2. Which uncertainty has the biggest effect on the estimated fluxes?

3. Derive an overall estimate for the accuracy of the Penman-Monteith method for both \(Q_E\) and \(Q_H\). Do this by assuming reasonable errors in the parameters and variables.

4. What is the impact of assuming neutral conditions versus actual stability conditions on the calculated value of \(r_a\)?

5. What error in the fluxes does this cause?

6. How will this error vary for extremely unstable or stable conditions?

Tasks: TODO

load necessary packages

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from pathlib import Path

load data

# code

examine data

# code and plot

Note

Please report issues or ask questions about this site on the GitHub page.

Data Sources

In this course, we focus on the AmeriFlux datasets for analysis. The URAO dataset might also be used, so some knowledge of both will be helpful.

Period selection

Note

Raw observational data undergoing a large amount of processing prior to be used (e.g. calculation of fluxes using Eddy Covariance techniques). They are usually collected with preliminary, or even without, “cleaning”.

It is thus crucial to consider different aspects for different parameters:

• cloud cover (e.g. this is important for the first plot)

• vegetation state (e.g. are leaves on the trees?) (LAI (Leaf Area Index))

• wind direction (e.g. does land cover type vary with wind direction?)

• Check that most of the data are available.

• Look carefully at the turbulent heat fluxes, i.e. small numbers of 30 min periods can be missing – that is ok

AmeriFlux

Tip

Become familiar with AmeriFlux website and data.

Sites for analysis

Group	Site Id	Name	Principal Investigator	Vegetation Abbreviation (IGBP)	Vegetation Description (IGBP)	Climate Class Abbreviation (Koeppen)	Climate Class Description (Koeppen)	Mean Average Precipitation (mm)	Mean Average Tempurature (degrees C)	Country	Latitude (degrees)	Longitude (degrees)	Elevation (m)	Years of Data
1	CA-Obs	Saskatchewan - Western Boreal, Mature Black Spruce	Andrew Black (andrew.black@ubc.ca)	ENF	Evergreen Needleleaf Forests	Dfc	Subarctic: severe winter, no dry season, cool summer	406	0.79	Canada	53.99	-105.12	628.94	1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010
1	US-Blk	Black Hills	Tilden Meyers (Tilden.Meyers@noaa.gov)	ENF	Evergreen Needleleaf Forests	Dfb	Warm Summer Continental: significant precipitation in all seasons	574	6.23	USA	44.16	-103.65	1718	2004, 2005, 2006, 2007, 2008
2	US-MMS	Morgan Monroe State Forest	Kim Novick (knovick@indiana.edu)	DBF	Deciduous Broadleaf Forests	Cfa	Humid Subtropical: mild with no dry season, hot summer	1032	10.85	USA	39.32	-86.41	275	1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018
2	US-MOz	Missouri Ozark Site	Jeffrey Wood (woodjd@missouri.edu)	DBF	Deciduous Broadleaf Forests	Cfa	Humid Subtropical: mild with no dry season, hot summer	986	12.11	USA	38.74	-92.2	219.4	2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015
3	US-Dia	Diablo	Sonia Wharton (wharton4@llnl.gov)	GRA	Grasslands	Csa	Mediterranean: mild with dry, hot summer	265	15.6	USA	37.68	-121.53	323	2010, 2011, 2012
3	US-KUT	KUOM Turfgrass Field	Joe McFadden (mcfadden@ucsb.edu)	GRA	Grasslands	Dfa	Humid Continental: humid with severe winter, no dry season, hot summer	777	7.9	USA	44.99	-93.19	301	2005, 2006, 2007, 2008, 2009
4	CA-Qcu	Quebec - Eastern Boreal, Black Spruce/Jack Pine Cutover	Hank A. Margolis (Hank.Margolis@sbf.ulaval.ca)	ENF	Evergreen Needleleaf Forests	Dfc	Subarctic: severe winter, no dry season, cool summer	950	0.13	Canada	49.27	-74.04	392.3	2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010
4	CA-Qfo	Quebec - Eastern Boreal, Mature Black Spruce	Hank A. Margolis (Hank.Margolis@sbf.ulaval.ca)	ENF	Evergreen Needleleaf Forests	Dfc	Subarctic: severe winter, no dry season, cool summer	962	-0.36	Canada	49.69	-74.34	382	2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010
5	US-Slt	Silas Little- New Jersey	Ken Clark (kennethclark@fs.fed.us)	DBF	Deciduous Broadleaf Forests	Dfa	Humid Continental: humid with severe winter, no dry season, hot summer	1138	11.04	USA	39.91	-74.6	30	2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014
5	US-UMB	Univ. of Mich. Biological Station	Christopher Gough (cmgough@vcu.edu)	DBF	Deciduous Broadleaf Forests	Dfb	Warm Summer Continental: significant precipitation in all seasons	803	5.83	USA	45.56	-84.71	234	2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018
6	US-Br3	Brooks Field Site 11- Ames	Tim Parkin (parkin@nstl.gov)	CRO	Croplands	Dfa	Humid Continental: humid with severe winter, no dry season, hot summer	847	8.9	USA	41.97	-93.69	313	2005, 2006, 2007, 2008, 2009, 2010, 2011
6	US-Bo1	Bondville	Tilden Meyers (Tilden.Meyers@noaa.gov)	CRO	Croplands	Dfa	Humid Continental: humid with severe winter, no dry season, hot summer	991	11.02	USA	40.01	-88.29	219	1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008
7	US-NC1	NC_Clearcut	Asko Noormets (noormets@tamu.edu)	OSH	Open Shrublands	Cfa	Humid Subtropical: mild with no dry season, hot summer	1320	16.6	USA	35.81	-76.71	5	2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012
7	US-Whs	Walnut Gulch Lucky Hills Shrub	Russ Scott (russ.scott@ars.usda.gov)	OSH	Open Shrublands	Bsk	Steppe: warm winter	320	17.6	USA	31.74	-110.05	1370	2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018
8	US-SRG	Santa Rita Grassland	Russell Scott (russ.scott@ars.usda.gov)	GRA	Grasslands	Bsk	Steppe: warm winter	420	17	USA	31.79	-110.83	1291	2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018
8	US-Var	Vaira Ranch- Ione	Dennis Baldocchi (Baldocchi@berkeley.edu)	GRA	Grasslands	Csa	Mediterranean: mild with dry, hot summer	559	15.8	USA	38.41	-120.95	129	2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018
9	US-Bo2	Bondville (companion site)	Carl Bernacchi (bernacch@uiuc.edu)	CRO	Croplands	Dfa	Humid Continental: humid with severe winter, no dry season, hot summer	991	11.02	USA	40.01	-88.29	219	2004, 2005, 2006, 2007, 2008
9	US-Br1	Brooks Field Site 10- Ames	Tim Parkin (parkin@nstl.gov)	CRO	Croplands	Dfa	Humid Continental: humid with severe winter, no dry season, hot summer	842	8.95	USA	41.97	-93.69	313	2005, 2006, 2007, 2008, 2009, 2010, 2011
10	CA-TPD	Ontario - Turkey Point Mature Deciduous	Altaf Arain (arainm@mcmaster.ca)	DBF	Deciduous Broadleaf Forests	Dfb	Warm Summer Continental: significant precipitation in all seasons	1036	8	Canada	42.64	-80.56	260	2012, 2013, 2014, 2015, 2016, 2017
10	US-Oho	Oak Openings	Jiquan Chen (jqchen@msu.edu)	DBF	Deciduous Broadleaf Forests	Dfa	Humid Continental: humid with severe winter, no dry season, hot summer	849	10.1	USA	41.55	-83.84	230	2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013

Tip

Refer to AmeriFlux Site Search page for full details of all sites used here.

Key variables for analysis

Tip

Refer to AmeriFlux Data Variables for full details of all measured variables.

• Surface Radiation Balance

Variable	Description	Unit
SW_IN	Shortwave radiation, incoming	\(\text{W m}^{-2}\)
SW_OUT	Shortwave radiation, outgoing	\(\text{W m}^{-2}\)
LW_IN	Longwave radiation, incoming	\(\text{W m}^{-2}\)
LW_OUT	Longwave radiation, outgoing	\(\text{W m}^{-2}\)

• Surface Energy Balance

Variable	Description	Unit
NETRAD	Net radiation	\(\text{W m}^{-2}\)
H	Sensible heat turbulent flux (no storage correction)	\(\text{W m}^{-2}\)
LE	Latent heat turbulent flux (no storage correction)	\(\text{W m}^{-2}\)
G	Soil heat flux	\(\text{W m}^{-2}\)

• Standard Meteorological Conditions

Variable	Description	Unit
TA	Air temperature	\(^{\circ} \text{C}\)
RH	Relative humidity, range 0-100	\(\text{%}\)
PA	Air pressure	\(\text{kPa}\)
P	Precipitation	\(\text{mm}\)
WS	Wind speed	\(\text{m s}^{-1}\)
WD	Wind direction	\({\circ}\)

• Others

Variable	Description	Unit
SWC	Soil water content (volumetric), range 0-100	\(\text{%}\)
LAI	Leaf area index by MODIS	\(\text{m}^{2} \text{ m}^{-2}\)

MODIS data associated with AmeriFlux sites

MODIS land products provides summaries of selected datasets for validation of models and remote-sensing products and to characterise field sites. Specifically, we use the LAI product in the class to help understand the relationships between surface parameters (e.g., albedo, surface resistance, etc.) and vegetation phenology.

Useful papers for understanding AmeriFlux datasets

1. Bergeron, Onil, Margolis, Hank A., Black, T. Andrew, Coursolle, Carole, Dunn, Allison L., Barr, Alan G., and Wofsy, Steven C. Comparison of carbon dioxide fluxes over three boreal black spruce forests in Canada. Global Change Biol, 13(1):89–107, January 2007. doi:10.1111/j.1365-2486.2006.01281.x.

2. Curtis, Peter S, Hanson, Paul J, Bolstad, Paul, Barford, Carol, Randolph, J.C, Schmid, H.P, and Wilson, Kell B. Biometric and eddy-covariance based estimates of annual carbon storage in five eastern north American deciduous forests. Agric. For. Meteorol., 113(1-4):3–19, December 2002. FLUXNET 2000 Synthesis. doi:10.1016/s0168-1923(02)00099-0.

3. Hollinger, D. Y., Ollinger, S. V., Richardson, A. D., Meyers, T. P., Dail, D. B., Martin, M. E., Scott, N. A., Arkebauer, T. J., Baldocchi, D. D., Clark, K. L., Curtis, P. S., Davis, K. J., Desai, A. R., Dragoni, D., Goulden, M. L., Gu, L., Katul, G. G., Pallardy, S. G., Paw U, K. T., Schmid, H. P., Stoy, P. C., Suyker, A. E., and Verma, S. B. Albedo estimates for land surface models and support for a new paradigm based on foliage nitrogen concentration. Global Change Biol., 16(2):696–710, February 2010. doi:10.1111/j.1365-2486.2009.02028.x.

4. Noormets, Asko, McNulty, Steve G., DeForest, Jared L., Sun, Ge, Li, Qinglin, and Chen, Jiquan. Drought during canopy development has lasting effect on annual carbon balance in a deciduous temperate forest. New Phytol., 179(3):818–828, August 2008. doi:10.1111/j.1469-8137.2008.02501.x.

URAO

Please visit URAO documentation site for information about observations at URAO.